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ABSTRACT 

In this paper, we combine the stellar spectral synthesis code STARBURST 99, 
the nebular modelling code MAPPINGS IHq, a 1-D dynamical evolution model of 
H II regions around massive clusters of young stars and a simplified model of syn- 
chrotron emissivity to produce purely theoretical self-consistent synthetic spec- 
tral energy distributions (SEDs) for (solar metallicity) starbursts lasting some 
10 8 years. These SEDs extend from the Lyman Limit to beyond 21 cm. We find 
that two ISM parameters control the form of the SED; the pressure in the diffuse 
phase of the ISM (or, equivalently, its density), and the molecular cloud dissipa- 
tion timescale. In particular, the shape of the FIR (dust re-emission) bump is 
strongly dependent on the mean pressure in the star-forming or starburst galaxy. 
This can explain the range of FIR colors seen in starburst galaxies. In the case 
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of objects of composite excitation, such diagrams potentially provide a means 
of estimating the fraction of the FIR emission that is contributed by an active 
nucleus. We present detailed SED fits to Arp 220 and NGC 6240, and we give 
the predicted colors for starburst galaxies derived from our models for the IRAS 
and the Spitzer Space Observatory MIPS and IRAC instruments. Our models 
reproduce the spread in observed colors of starburst galaxies. From both the 
SED fits and the colorxolor diagrams, we infer the presence of a population of 
compact and ultra-compact H II regions around single OB stars or small OB 
clusters. Finally, we present absolute calibrations to convert observed fluxes into 
star formation rates in the UV (GALEX), at optical wavelengths (Ha), and in 
the IR (IRAS or the Spitzer Space Observatory). We show that 25/im fluxes are 
particularly valuable as star formation indicators since they largely eliminate one 
of the parameters controlling the IR SED. 

Subject headings: ISM: dust — extinction, H II-galaxies:general,star formation 
rates, starburst-infrared:galaxies — radio continuum:galaxies — ultraviolet:galaxies 

1. Introduction 

A knowledge of the star formation rate (SFR) is fundamental if we are to understand 
the formation and evolution of galaxies. In a much-cited seminal paper Madau et al. (1996) 
connected the star formation in the distant universe to that estimated from low-redshift 
surveys, by plotting the estimated star formation rate per unit co-moving volume against 
redshift. A wide variety of techniques contribute towards populating this diagram with 
observational points. 

At the present time an unacceptable degree of uncertainty still attaches to the Madou 
plot, and hence to our overall understanding of the evolution of star formation in the universe. 
This is largely the result of the effects of dust obscuration in and around the star-forming 
regions, which are particularly severe for those techniques which depend upon optical or UV 
data (see the review by Calzetti 2001). Use of the far-IR may provide improved estimates 
(Dwek 1998; Gispert et al. 2000). Nonetheless, uncertainties still remain, since to derive 
these rates, one must assume that the dust is effectively acting as a bolometer wrapped 
around the star forming region, and that cool and old stars do not provide too much of a 
contribution. All of these assumptions deserve closer examination in both disk and starburst 
galaxies. 

A further major uncertainty relates to the fact that we do not yet have reliable and 
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self-consistent theoretical spectral energy distributions (SEDs) covering the full range of the 
electromagnetic spectrum. In principle, almost any part of the SED of a starburst can be 
used as a star formation indicator, provided that the appropriate bolometric correction to 
the absolute luminosity can be made. These bolometric corrections critically depend on the 
dust absorption, and the geometry of the dust clouds, both optically thin and optically thick, 
with respect to both the ionizing stars and the older stellar population. 

In principle, provided that the initial mass function (IMF) is invariant, the intrinsic UV 
luminosity should scale directly as the star formation rate. From stellar spectral synthesis 
models, the star formation rate is given in terms of the 1500A flux by (Kennicutt 1998; 
Panuzzo et al. 2003): 
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In practice, the results obtained with this technique depend critically on the corrections 
needed to account for the dust in and around the star-forming galaxy. These are very 
uncertain; there are claims that a typical z = 3 galaxy suffers a factor of 10 extinction at 
a wavelength of 1500A in the rest -frame of the galaxy (Meurer, Heckman & Calzetti 1999; 
Sawaki & Yee 1998). Others argue for more modest corrections (e.g. Trager et al. 1997). In 
normal galaxies, Bell & Kennicutt (2001) find that although mean extinctions are modest 
(about 1.4 mag at 1550A), the scatter from galaxy to galaxy is very large. For individual 
galaxies at high redshift the reddening and the SFR are strongly correlated, and since many 
fainter galaxies at high redshift are missed completely, potentially large uncertainties exist 
in the derived star formation rates (Adelberger & Steidel 2000). 

A fairly direct and extensively-used technique is to measure hydrogen recombination line 
fluxes (Bell & Kennicutt 2001; Gallego et al. 1995; Jones & Bland-Hawthorn 2001; Moorwood 
et al. 2000; Tresse & Maddox 1998; Pascual et al. 2001; Yan et al. 1999; Dopita et al. 2002, 
to name but a few). Provided that the H II region can absorb all the EUV photons produced 
by the central stars, this should be a reliable technique, since in this case the flux in any 
hydrogen line is simply proportional to the number of photons produced by the hot stars, 
which is in turn proportional to the birthrate of massive stars. This relationship has been 
well calibrated at solar metallicity for the Ha line. In units of M yr -1 , the estimated star 
formation rate is given by (Dopita & Ryder 1994; Kennicutt 1998; Panuzzo et al. 2003): 
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Provided that the star-formation regions are resolved, their Balmer decrements can be 
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used to estimate the absorption in any foreground dust screen. However, it is possible that 
some star forming regions are completely obscured, even at Ha. This problem can be avoided 
by observing at infrared wavelengths, in Bra, for example. A further complication is that the 
dust content of the nebula is metallicity-dependent, and therefore an appreciable fraction of 
the hydrogen ionizing photons may be absorbed by dust in high metallicity H II regions. This 
possibility, first seriously quantified by Petrosian, Silk & Field (1972), has been discussed 
by a number of authors since (Panagia 1974; Mezger, Smith, & Churchwell 1974; Natta & 
Panagia 1976; Sarazin 1977; Smith, Biermann, & Mezger 1978; Shields & Kennicutt 1995; 
Bottorff et al. 1998). Its effect has been investigated and quantified (as far as is possible 
by direct observation) in a series of recent papers by Inoue and his collaborators (Inoue, 
Hirashita, & Kamaya 2000; Inoue, Hirashita & Kamaya 2001; Inoue 2001). Dopita et al. 
(2003) has shown that this effect increases in importance in the more compact H II regions, 
and may lead to errors as high as ~ 30% in the estimated SFR. 

In the case of dusty starburst galaxies, the total 8-1000 /im infrared luminosity, Lir, as 
measured in the rest frame frequency of the galaxy can be used to estimate the total SFR 
(Kennicutt 1998; Panuzzo et al. 2003): 

-^t1 ■ (3) 
erg.s 1 

In deriving the proportionality between IR flux and SFR, it is generally assumed that 
the dust acts as a sort of "bolometer", wrapped around the star formation regions, and 
re-emitting most, if not all, the incident EUV and FUV photons. In practice, this is one 
of the many "geometrical" factors which render the derivation of the SFR quite uncertain. 
For example, Inoue (2002) includes a number of these geometrical factors in the following 
formula for the IR luminosity of a galaxy: 
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Lm = L Lya + (1 - f)L LC + eL uv + r]L old (4) 

This allows a fraction, (1 — /), of the ionizing flux to be absorbed by dust in the 
H II region, assumes that all of the Lymana photons are multiply scattered by the gas and 
ultimately absorbed by dust in the surrounding H I region, that a fraction, e, of the non- 
ionizing UV photons are absorbed, and that a fraction rj of the radiation field produced by 
the old stars is also absorbed by dust. Photoionization models of H II regions show that, 
typically, the Lymana flux is of order 30% of the total stellar flux absorbed by the nebula, 
so that globally, the dust is fairly efficient in capturing and re-radiating the flux originally 
radiated by the star in the Lyman continuum. Indeed, since most of the FIR flux results from 
the absorption of both the non-ionizing and ionizing photons in the photodissociation regions 
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around H II regions, and since the temperature of the dust depends critically on the distance 
of the dust from the stars, it is clear that a clear understanding of the geometry of the dusty 
H I and molecular gas with respect to newly-born stars is critical to our understanding of 
the SED. 

In combination with STARBURST 99 (Leitherer et al. 1999) modelling, Hirashita, Buat 
& Inoue (2003) have applied equation 4 to the observations of star forming galaxies in order 
to estimate the empirical values for the geometrical factors, which provide estimates for SFRs 
consistent between the UV, IR and Ha estimators. 

Recently, several very sophisticated approaches to the determination of the SED in 
galaxies have been developed. For normal galaxies the models of Silva et al. (1998) explicitly 
allow for the clumpy distribution of dust, strongly spatially correlated with the UV-emitting 
young stellar clusters, and probably associated with the dense molecular clouds from which 
the young stars have been formed. The model of Popescu et al. (2000), subsequently extended 
in later papers (Misiriotis et al. 2001; Tuffs et al. 2004), explicitly takes into account the 
three-dimensional stucture of spiral galaxies both in terms of their stellar and gas and dust 
distributions to provide projected photometric properties and attenuation characteristics. 

For starburst galaxies, Hirashita & Ferrara (2002) and Takeuchi et al. (2003) have 
developed a dust re-emission model appropriate to low-metallicity young galaxies at high 
redshift. Tagaki, Vansevicius & Arimoto (2003) and Tagaki, Arimoto & Hanami (2003) have 
developed a mathematically elegant spherically symmetric radiative transfer model which 
includes the very important temperature fluctuations of the small grain populations and 
emission in the PAH features between 3 and 30/im. Applied to starburst systems, this model 
suggests that star formation rates are positively correlated with the structural characteristics 
of the star forming region: regions of high star formation rate are both compact and highly 
obscured at optical wavelengths. 

Because of their effect on both the efficiency of the production of IR emission and on 
the form of the IR emission through the dust temperature distribution, these geometrical 
factors largely determine the far-IR SED. In a pair of recent papers Dale and his co-workers 
have used ISO observations to develop a (largely empirical) family of templates to fit the 
variety of the forms of the IR SED (Dale et al. 2001; Dale & Helou 2002). This work assumes 
a power law distribution of dust mass with local radiation field intensity to provide a wide 
range of dust temperatures. It appears to indicate that the SEDs of star forming galaxies 
could fit to a one-parameter family of curves determined essentially by the exponent, a, of 
the power law distribution of the dust masses with the radiation field intensity, where the 
radiation field is assumed to have the color of the local ISRF. This will not be the case in 
starburst galaxies, where the radiation field will be harder than the local ISRF, and where 
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the geometry of the dust with respect to the gas is determined by the evolution of H II 
regions. Also, it may not be applicable even in normal galaxies, which are known through 
direct observation to appear as the superposition of different geometrical components whose 
dust emission has different intensities and colors (Tuffs & Popescu 2003). The commonly 
used 60/100 /im color diagnostic, cannot be interpreted in terms of a single grain emission 
component, since the 60/xm emission is seen to come mainly from warm dust locally heated 
in HII regions, the 100/xm emission comes both from the HII regions and the cold diffuse 
inter-arm emission. 

In this paper, we attempt to provide a theoretical approach to the generation of pan- 
spectral energy distributions in the context of starburst galaxies. This provides physical 
insight into how the geometrical factors discussed above are determined largely by the tem- 
poral evolution of the H II region complexes in starburst galaxies. We have included the dust 
physics and radiative transfer both in the ionized H II plasma, and in the dense molecular 
clouds associated with them. In particular, we will show that the stochastic average pressure 
within the star formation region plays a major role in determining the FIR SED through 
its effect on the grain temperature, and that the rate at which the molecular clouds are dis- 
persed and destroyed by the expanding H II complexes is fundamental for our understanding 
of the geometrical factors which affect the global SED. 

This paper is organized as follows. In Section 2 we describe the STARBURST 99 
modelling which provides the input stellar spectra, and in Section 3 we uses the mechanical 
energy input given by the STARBURST 99 modelling to build a simple one-dimensional 
temporal evolution model for the H II regions surrounding a typical OB stellar cluster. In 
Section 4 we describe the dust model used in this model, and show how the photodissociation 
of PAH molecules combined with a turbulent ISM may explain the observed attenuation curve 
of starburst galaxies. The MAPPINGS Illq (Sutherland & Dopita 1993; Dopita et al. 2002) 
modelling of the emission line, atomic continuum, and dust re-emission spectrum of these 
H II regions and their surrounding photodissociation regions (PDRs) and the generation 
from these models of the global starburst SEDs is described in Section 7. In Section 8 we 
move on to discuss the parameters which have a controlling influence of the form of the 
theoretical SEDs, and go on to compare the models with observations in Section 9. Finally, 
in Section 10 we use the model to provide a theoretical re- calibration of the various star 
formation indicators. In the conclusion, we enumerate the main results of the paper, and we 
emphasize the limitations of the models. 

All the models presented in this paper have a solar abundance set. In Paper II of this 
series, we will investigate the effect of chemical abundance on the form of the SEDs. 
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2. Stellar Spectral Synthesis Modelling 

Starburst activity was first linked to galaxy interactions decades ago. In 1967, Sersic & 
Pastoriza (1967) found that galaxies harbouring hotspots all had extreme blue colours in their 
centers, a property shared by the majority of interacting or peculiar galaxies in the Vorontsov- 
Velyaminov (1959) and Arp catalogues. An excess of radio emission was found in interacting 
galaxies, strengthening further the link between star formation and interaction (Sulentic 
1976; Condon & Dressel 1978; Hummel 1980). Larson & Tinsley (1978) first suggested 
that tidal forces in interacting galaxies could trigger bursts of star formation. Since then, 
numerous studies have produced evidence for interaction-induced starburst activity (eg., 
Petrosian 1978; Condon et al. 1982; Keel et al. 1985; Madore 1986; Kennicutt et al. 1987; 
Bushouse 1987; Hummel et al. 1990; Donzelli & Pastoriza 1997; Barton, Geller, & Kenyon 
2000; Petrosian et al. 2002; Lambas et al. 2003). 

The typical duration of such interactions is a dynamical timescale; comparable to the 
orbital timescale of stars in the gravitational well of the system, ~ 10 8 years. The details 
of the star formation triggered by such interactions depend on the many parameters of the 
interaction such as the impact parameter, the relative velocities, the relative orientations of 
the spin axes, and the relative masses of the systems. Simulation of galactic mergers (Barnes 
& Hernquist 1996) suggests that the most violent star formation tends to be triggered in 
two separate episodes; during the initial collision of the systems and later, when tidal tail 
debris rains down into the nuclear regions of the merged systems. Enhanced star formation 
can occur throughout the interaction. 

Given that each case will be different, for simplicity in what follows, we have assumed 
that in our model starbursts, stars are being created continuously as separated 10 4 M 
clusters scattered through the same volume, so as to give a mean star formation rate of 
1 M yr~ x . Thus, in the total lifetime of the starburst, 10 s years, some 10 4 such clusters will 
be formed. 

Clusters of young stars are dust-enshrouded for the first ~ 1 — 2 Myr, and produce the 
majority of their ionizing UV photons within the first ~ 4 — 6 Myr of their lifetime. Thus 
star formation indicators based upon nebular emission lines such as Ha or [O II], or else 
based upon dust re-emission processes such as the 60 /xm continuum flux are strongly biased 
towards measuring recent star formation rates. By contrast, features which more directly 
measure the older population or the total luminosity, such as the near-IR continuum or the 
PAH emission features measure the star formation rates integrated over somewhat longer 
periods. Assumptions about the total duration of the starburst do not greatly affect the 
calibration of the star formation rates based on either the UV, emission line, or FIR SED - 
the total luminosity of starburst increases by only by 10% between 5 x 10 7 years and 10 s years 
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according to the STARBURST 99 models (Leitherer et al. 1999). 

We have used the most recent version of the STARBURST 99 code to generate the 
theoretical stellar SEDs of: 

• Instantaneous starbursts, sampled at 1 Myr age intervals from Myr up to 10 Myr 
and 

• Continuous star formation at the rate of 1 M yr _1 up to a maximum age of 10 8 years. 

These models assume solar metallicity, a standard Salpeter initial mass function, and 
upper and lower mass cutoffs of 120 M and 1 M Q respectively on the zero-age main sequence. 
The stellar atmosphere models for stars with plane parallel atmospheres are based on the 
Kurucz (1992) models as compiled by Lejeune, Cuisinier & Buser (1997). The fully line- 
blanketed Wolf-Rayet atmosphere models of Hillier & Miller (1998) and the non-LTE O- 
star atmospheres of Pauldrach, Hoffman, & Lennon (2001) have been incorporated into 
STARBURST 99 as described in Smith, Norris, & Crowther (2002). We use the full suite 
of atmosphere models to provide the most up-to-date stellar radiation field possible. All 
STARBURST 99 models were re-binned to the lower spectral resolution input energy bins 
required by MAPPINGS Illq. 

The instantaneous burst models were used as the input spectra for the H II region neb- 
ular gas and dust modelling, while the continuous star formation models form a theoretical 
global SED for the older stars with ages between 11 and 100 Myr. The contribution of the 
older stars was distributed between the H II regions of different ages in proportion to their 
volume, and added to the stellar SEDs of the young cluster which principally excite the H II 
regions. Thus our global SED is composed of the sum of all the contributions from the stars 
and their surrounding H II regions which provide an average star formation rate of 1 M yr _1 
for a total time of 10 8 years, finely sampled in time at 1 Myr time intervals between and 
10 Myr. 

3. The Dynamical Evolution of H II Regions 

It is very easy to understand why the temporal evolution of the H II regions surrounding 
young clusters of stars might be quite important in influencing the SED of starburst galaxies. 
First, the inner radius, i? in , of the ionized shell and the pressure in the H II region together 
determine the ionization parameter, U = F^/AncR^nmi, where F* is the flux of ionizing 
photons from cluster stars, c is the speed of light and nnn is the particle density in the 
H II region. As is well understood from photoionization theory, see Dopita & Sutherland 
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(2003), the emission line spectrum of an H II region is determined primarily by the effective 
temperature of the cluster stars (itself a function of time) and by the ionization parameter, 
which is determined by the radius : time, and the radius : pressure relationships. These 
relationships need to be derived from models of the dynamical evolution of the H II regions, 
and are determined by the deposition of mechanical energy and momentum from the stellar 
winds and by supernovae. 

The FIR SED is also influenced by these parameters, because, provided the non-ionizing 
UV is absorbed on a scale less than the radius of the HII region, the outer radius of the H II 
region largely determines the characteristic dust temperature for a given stellar luminosity, 
and the IR PAH features are generated mostly in the photodissociation region surrounding 
the H II region. 

In essence, then, the evolution of the H II regions determines the geometry of the gas and 
the dust with respect to the hot stars, and these geometrical effects determine various aspects 
of the SED. These effects have not hitherto been taken into account in any self-consistent 
way in the computation of SEDs. The most sophisticated treatment so far can be found in 
the "semi-analytic" framework of Granato et al. (2000); Bressan, Silva & Granato (2002) 
and Panuzzo et al. (2003). In this, the young stars are treated as escaping from a dense 
molecular cloud environment on a characteristic time scale which is environment dependent. 
However, such a treatment cannot properly take into account either the change of ionization 
parameter within the H II regions, nor the change in the mean dust temperature with time. 

The theory of the size distribution of H II regions has been presented by Oey & Clarke 
(1997) and Oey & Clarke (1998). According to these papers, the H II regions expand and 
evolve as mass-loss and supernova driven bubbles for as long as their internal pressure exceeds 
the ambient pressure in the ISM. When this condition is no longer met, the H II region is 
assumed to "stall". Of course, on reaching the stall condition, the H II region does not 
abruptly stop expanding, but a momentum-conserving expansion continues until interstellar 
turbulence destroys the integrity of swept-up shell. The time taken to reach the stall radius 
turns out to be proportional to the stall radius, and the mass of the OB star cluster is 
important in determining the stall radius; low mass cluster H II regions stall early and at 
small radius. 

This theory emphasizes the important role of the pressure in the ISM in determining the 
radius of an individual H II region, and its evolution with time. In particular, with higher ISM 
pressure, the H II region stalls at smaller radius, and therefore the mean dust temperature 
in the atomic and molecular shell around the H II region should be higher in high-pressure 
environments. We would therefore expect the FIR SED to peak at shorter wavelengths in 
high-pressure starburst regions. It is this correlation that we seek to investigate in this paper. 
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According to the mass-loss driven bubble theory of Weaver et al. (1977), the radius, R, 
and the internal pressure, P, of an H II region are given by (c.f. equations 24 and 25 of Oey 
& Clarke (1997)); 

fi=(j^y /5 ^» i/5 < 3/5 <5) 
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Here, L mcc h is the mechanical luminosity of the mass-loss from the central stars, po is the 
density of the ambient medium and t is the time. The particle density is given in terms of 
the ambient density by n = po/^rriH, and the ambient pressure Po — nkTy. Eliminating t 
between equation 5 and 6, and setting P = P we obtain the following conditions for stall: 
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These equations imply that, for any given pressure in the ISM, all stalled H II regions 
have a common ratio L mcc h/nHii-R 2 . Thus to the extent that oc L me ch, where F* is the flux 
of ionizing photons from the central stars, then all stalled H II regions will be characterized by 
a common value of the ionization parameter, U = F # /4ncR?jiiw- 111 fact, the STARBURST 
99 models show that the ratio of mechanical luminosity to ionizing photons may vary by as 
much as a factor of ten during the first 10 7 years of evolution. Nonetheless, most observed 
extragalactic H II regions show a fairly limited range in their ionization parameters (Kewley 
& Dopita 2002), and most H II regions should be in the stall condition (Oey & Clarke 1997), 
so that the above argument may well represent the basis of a theoretical explanation of this 
curious observational result. 

Although we have emphasized the role of pressure, equation 8 shows that the stall 
condition might be better cast as a relationship between radius and ambient density, since 
it is the ambient density that determines the expansion history of the bubble. The effective 
ambient density is set by the density of the phase of the ISM which has the largest volume 
filling factor. We have assumed that, in starburst galaxies, this phase is a warm ionized 
medium (WIM) or warm neutral medium (WNM) heated by the photons and mechanical 
energy of the starburst itself. Since both of these phases have a temperature T ~ 10 4 K, they 
also have similar densities for a given pressure. 



Whether or not the H II regions are in the stall condition, the observed ionization 
parameter provides an important constraint on the internal pressure, and, more importantly, 
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on the radiation field density incident on photodissociation region (PDR) around the H II 
region. This determines the characteristic dust temperature, and hence the wavelength of 
the peak in the FIR dust re-emission feature. 

H II regions in the stalled condition are found preferentially in regions of higher ISM 
pressure or density, and around lower mass clusters. Here, we consider the variation in 
pressure, but we have not taken into account the mass distribution of the central clusters, 
using instead a 'representative' cluster of mass 1O 4 M . Clearly, the cluster mass distribution 
is an important parameter, the investigation of which is deferred for a future paper. 

In the stall condition the expansion of the bubble does not cease entirely. By the 
time the stall radius is reached there is already a massive shell of swept-up gas, and this 
carries a substantial momentum. In our evolutionary calculations, we have assumed that the 
evolution of the bubble subsequent to the stall condition being reached continues according 
to the momentum conserving solution for the swept-up shell: 



R = 
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where the mass and velocity of the swept up shell at the stall radius are M sta ii and i> s taii, 
respectively. 



3.1. Bubbles with Time Dependent Luminosity 

We have used the stellar wind and supernova power output for a 1O 4 M cluster as 
predicted by the STARBURST 99 code to make a Runge-Kutta integration of the equation 
of motion of the H II shell. The internal pressure is calculated from equation 8, with the 
instantaneous power output being smoothed over a period of lMyr to eliminate short time- 
scale pressure fluctuations. 

For a time dependent wind + supernova mechanical luminosity L m cch(£) and an ambient 
medium density of po, the equations of conservation of mass, momentum and energy for the 
mass M of the swept-up shell for a stellar wind bubble are; 



j f [M] = A7TR 2 p R, (10) 

4 [MR] = 4nPR 2 , (11) 
at 
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where R is the outer shell radius. Here P is the bubble pressure, and P » Pq. By assuming 
that the bubble pressure is much greater than the ambient pressure, we have an energy driven 
bubble. This will be true until shortly before the stall radius is reached. 

Assuming that 7 = 5/3 defines the equation of state of the bubble gas, we can combine 
these into a single differential equation, which is numerically integrated to allow for an 
arbitrary, time dependent, mechanical luminosity L w , 



This equation was then integrated using a standard Runge-Kutta 4th order integration with 
an adaptive step size, using the STARBURST 99 generated tabulated input mechanical 
energy luminosity for our cluster as a function of time until the stall condition was reached, 
after which the momentum-conserving solution was adopted (equation 9). 

This simple integration does not yield results which are in accord with observations. 
The bubble expands to too large a diameter without reaching the stall condition, and the 
internal pressure is too high. As a result, the ionization parameter in the ionized shell is 
about a factor of ten below the observed values for extragalactic H II regions as determined 
by Kewley & Dopita (2002) at all times during the lifetime of the bubble. 

This problem can be rectified by assuming that the mechanical power output of the 
young stars is much less effective than the simple 1 — D hydrodynamical models would 
suggest in terms of its efficiency in inflating the ionized bubble. Empirically, we determined 
that the effective power input of the stars would have to be decreased by a factor of ten 
to ensure that the ionization parameter in the H II region matches the observed range of 
3 x 10 -4 < U < 3 x 10~ 3 (Kewley & Dopita 2002). This is in part due to the decrease in the 
radius of the bubble, and in part due to the lower internal pressure. 

Interestingly enough, similar discrepancies between observation an the simple mass-loss 
bubble theory have emerged in analyses of the dynamics of bubbles and of their stellar con- 
tent ( Garcia- Segura et al. 1995; Oey 1996; Naze et al. 2001). In these studies discrepancies 
of between one and two orders of magnitude have been noted. In all cases, the energy flux 
in the stellar winds predicted by the analysis of the stellar content exceeds the energy flux 
inferred from the dynamics. These discrepancies probably arise out of the limitations of the 
1 — D hydrodynamic models. Such models cannot take into account a number of effects 
which lead to a decreased efficiency in inflating the mass-loss bubble. These include en- 
hanced radiative losses in the shocked stellar wind plasma through operation of the Vishniac 
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instability, blowouts caused by the breakup of the atomic and molecular shell surrounding 
the H II region, and 'poisoning' of the hot plasma through photo-ablation and entrainment 
of fragments within the hot shell. 

In the light of this, we have modelled our 1O 4 M clusters as if they only had the 
mechanical energy luminosity of a 10 3 M Q cluster. This ensures that our model H II regions 
have ionization parameters which accord with observations, and this will also ensure that 
the dust grains in the PDRs of these bubbles will be illuminated by a radiation field of the 
correct strength, so that the grains will have the correct temperatures. 

In figure 1 we plot the computed time evolution of the radius of the bubble as a function 
of the mean pressure P/k (measured in units of cm~ 3 K). The undisturbed gas outside the 
bubble is assumed to be part of the warm ionized medium, so that P/k = 10 4 corresponds 
to an atomic density of ~ 1.0 cm" 3 . In figure 2 we plot the corresponding evolution of the 
internal pressure in the bubble, which can be identified as the pressure in the H II region. 
Note the peak in the pressure which appears near log(t)~ 6.5 yr, which corresponds to the 
onset of the energy input by Type II supernovae. 

For P/k < 10 6 the bubble does not stall until quite late in its evolution, but for the 
highest pressures the stall condition is reached very early on. This ensures a large range 
in the final radii. Furthermore, for all the models,the pressure inside the bubble is quite 
strongly correlated with the external pressure, for times greater than about lMyr and up 
to the time at which the Type II supernova energy input switches off. This suggests that a 
measurement of the density in the ionized plasma using a density-sensitive ratio such as the 
[S II] AA6717/673lA ratio could be used to help estimate the mean ambient ISM pressure 
in the star-forming region. 

During the evolution of the H II bubble, the distribution of the radiation field within 
the bubble can be modified by the dynamical evolution of the central cluster, and indeed, 
by the turbulent break-up of the shell of gas around the H II region. For, example, at 10 7 yr 
the mean expansion speed of the H II bubble is 0.6, 5.7 and 17.3 kms" 1 for Po/k = 10 8 , 10 6 
and 10 4 cm _3 K, respectively. Thus, at the higher pressures, the velocity dispersion of the 
stars in the central star cluster is likely to rival or exceed the expansion of the H II region, 
leading to a breakdown in the assumption that the illumination of the PDR at the edge of 
the bubble is from a single source at the center of the bubble. Under such circumstances, 
the rate of expansion of the bubble may be determined by the velocity dispersion of the stars 

Fig. 1. — The time variation of the radius of the H II bubble as a function of the pressure in 
the ISM 
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rather than the hydrodynamical effects, with the "bubble" in fact being a cluster of smaller 
H II regions from individual stars, rather than a single cavity. In this case, the picture might 
be of a molecular cloud being excavated from within through the growth of independent HII 
regions inflated by single stars (or by small numbers of stars), eventually creating a porous 
parent molecular cloud. 

All of this serves to highlight the fact that the time has come to abandon the simple 
1-D evolutionary hydrodynamic model in favour of an approach which properly takes into 
account both the dynamics of the central cluster and the hydrodynamics of both the ionized 
and the molecular gas. This treatment is beyond the scope of this paper. 

4. Treatment of Dust and PAH Physics 

4.1. The PAH Molecules 

The numerous emission features in the 3 — 30/im region, formerly referred to as the 
Unidentified Infrared Emission Bands or UIBs were first identified as due to polycyclic aro- 
matic hydrocarbons (PAHs) by Leger & Puget (1984). Since that time, our understanding of 
these features has become rather sophisticated, as more and more observational constraints 
have been produced, principally by the ISO mission (Peeters 2002). At the time of writing, 
the Spitzer Space Observatory is providing stunning new images in the PAH bands. Clearly, 
any serious attempt to model SEDs must now include a good deal of PAH physics, and in 
the following subsections we describe the treatment of PAHs within the mappings code. 

In modelling the PAHs within photoionized regions we have opted for simplicity over 
an exact treatment, especially as the constraints on both the physics and the species of 
PAHs remain so poor. In the MAPPINGS Illq code we have simply assumed that all the 
PAHs can be represented by a single representative type, coronene (C24H 12 ). Not only has 
coronene been observed within the ISM, but it is also pericondensed, reasonably stable and 
can adequately match the observational constraints in the UV and IR. Though perhaps 
slightly smaller than the average ISM PAH (Allamandola, Hudgins, & Sandford 1999), it 
has the great advantage that there exists a large amount of laboratory data for it, hence 
any modelling done with these data can be considered reasonably accurate from a purely 
physical viewpoint. 

Fig. 2. — The time variation of the interior pressure of the H II bubble as a function of the 
pressure in the ISM 
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The spectroscopic properties (emission and absorption) of a PAH molecule depends 
upon the number of carbon atoms (Nc), the hydrogenation (H per C) and the charge state 
of the molecule. In addition to these, the optical properties are affected by the shape of 
the PAH (cata- or peri- condensed) and the addition of end groups (like methyl -CH 3 ). In 
the modelling, we assume a pericondensed structure with no additional end groups. The 
PAH abundance within the code is determined by the amount of Carbon depleted onto PAH 
molecules. This was chosen to provide a good match to the absolute strength of the PAH 
emission features in the infrared. 

The use of a single species of PAH makes the determination of PAH survival, optical 
properties and charge much easier to determine. The absorption, emission and photoelectric 
processes on a single size and shape are also much more readily calculated. Li & Draine 
(2001), using all this available data have constructed absorption cross sections for both 
neutral and ionized "astronomical" PAHs from the FUV to FIR. They did this for a range 
in N c , with the largest of PAHs blended in with the optical properties of graphite grains. 

Within MAPPINGS Illq, we have used their data 1 to represent the PAH optical prop- 
erties. The PAH optical data of Li & Draine (2001) is defined in terms of a PAH "radius", 
such that there is a continuum in data between PAH and graphite grains. This PAH radius 
connects directly with Nc, with a = 10(Nc/468) 1 / 3 A. This gives for coronene an effective 
radius of a = 3.7A and surface area of 173A 2 . These values, and the corresponding Li 
& Draine (2001) data for this size, are used within the code to define the effective PAH 
molecule. We use both the neutral and ionized PAH optical data for the code. The fraction 
contributed to the total PAH opacity is determined by the charge state of the model PAH, 
which is determined by the balance of electron sticking and photoelectric ejection. 



4.2. The PAH Emission Spectrum 

The emission spectrum of PAHs is determined by the natural modes of vibration, bend- 
ing and other deformation of the carbon skeleton, modified by the effect of the dangling 
groups and the electric charge state of the molecule. The principal modes and their charac- 
teristic wavelengths are as follows: 

• 3.29/iin Aromatic C-H stretch (v — 1 — 0); 

• 3.40/im Aliphatic C-H stretch (the v — 1 — anti-symmetric mode); 



Available at http://www. astro. princeton.edu/^draine /dust/dust. did. html 
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• 3.51/j.m Aliphatic C-H stretch (the v — 1 — symmetric mode); 

• 6.2yum and 7.7-7.9/iinC-C skeletal deformations; 

• 8.6/iin C-H in-plane bend; 

• 11.2/im C-H out-of-plane bend (solo mode); 

• 11.9//m C-H out-of-plane bend (duo mode); 

• 12.7/im C-H out-of-plane bend (tri mode); 

• 15-20/im C-C-C bending modes (Van Kerckhoven et al. 2000). 

These emission features can be adequately represented by Lorentzian fits (Boulanger et al. 
1998) of the form: 

F(x) = - r {A/7Ta \ n (14) 

[l + (x- x ) 2 /a*} 

in which x — 1/A cm -1 , the central wavenumber of the feature is Xq, the FWHM = 2cr, and 
the area under the curve is A. 

Two model fits were made, the first to the data of Boulanger et al. (1998) for the Rho 
Ophiuchus dark cloud and for NGC7023, and the second to the data of Laurent et al. (2000) 
for the M82 starburst. The second fit used more components, as listed in Verstraete al. 
(2001), and also included a 2 component fit to the 15-20/im C-C-C bending modes (Van 
Kerckhoven et al. 2000). This fit applies to the Orion Bar region, M17-SW, and NGC2023. 
The final model PAH emission spectrum adopted for the MAPPINGS Illq code is given in 
Table 1. 

Emission within the MAPPINGS Illq code is treated as an energy conservation process. 
In equilibrium all the energy gained by a PAH through the absorption of photons can be either 
lost through photoelectric processes or IR emission. The fraction of the energy lost through 
photoelectric processes is determined, and we assume the remaining fraction of energy is 
re-emitted in the IR according to our empirical fit to the astrophysical PAH emission bands 
(UIBs). 

This is not an exact treatment as the PAH molecules will likely undergo stochastic 
heating processes and will lose some of their energy through IR continuum emission instead 
of via these fluorescent bands. However the code also allows for very small graphite grains 
which do undergo stochastic heating, thus both forms of emission are catered for. 



Table 1. Fit to PAH Emission Components 



Parameters: 



A (/im) 


3.39 


6.2 (#1) 


6.2 (#2) 


6.2 (#3) 




on/in o 

3040.3 


1608.5 


1608.5 


1593.9 


Peak 


1.049E-4 


a r—r(~\ r~\ T — \ A 

4.738E-4 


6.768E-4 


3.113E-4 


a 


20.40 


31.62 


12.1 


34.90 


A (/iin) 


6.949 


7.5 


7.6 


7.8 


x 


1 A A C\ 

1440 


i o o o r - 

1328.5 


1311.3 


1274.3 


Peak 


i A , "\ 1 A T — I A 

3.214E-4 


6.768E-4 


9.730E-4 


1.438E-3 


a 


100 


44.72 


14.0 


35.0 


A (/iin) 


8.6 


10.0 


11.3 (#1) 


11.3 (#2) 


x 


1163.1 


998 


889.4 


889.8 


Peak 


1.015E-3 


2.115E-4 


1.015E-3 


2.200E-3 


a 


30.0 


158 


8.8 


5.68 


A (/im) 


H.3 (#3) 


12.7 


13.6 


14.3 


x 


880.9 


787 


735 


699 


Peak 


1.354E-3 


1.015E-3 


1.184E-4 


1.099E-4 


a 


15 


30.1 


14.9 


14.9 


A (/im) 


16.4 


17.8 






x 


610 


562 






Peak 


4.230E-4 


2.961E-4 






a 


17 


30.0 
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4.3. PAH Photodestruction 

Where detailed measurements have been made of the distribution of PAH emission 
within compact H II regions, it is found that the PAH avoids the ionized regions, and is 
instead found in a narrow zone of emission (the photo-dissociation region) beyond the outer 
boundary of the ionized gas (Burton et al. 2000). This is clear evidence that PAH molecules 
cannot survive for significant lengths of time in the hostile environment in the ionized zones 
of an H II region. This is to be expected, because the photodissociation timescales in this 
region are very short. However, PAHs can survive, and are excited into emission in the 
photodissociation regions (PDRs) outside the ionized H II region. 

The condition for photodestruction may be set in a number of ways. The simplest is 
that the presence of hard photons leads to rapid photodissociation in the ionized gas. In this 
case, the PAHs are destroyed in the ionized gas and their C is returned to the gas phase. 

A somewhat more sophisticated approach would be to say that PAHs are destroyed 
whenever the absorption weighted mean radiation field (not necessarily ionizing) results in 
a photodissociation timescale which is short compared to the residence time of the PAHs 
in that radiation field. In this case we would compare photodissociation timescales with 
dynamical timescales to obtain a criterion for survival. 

A still more sophisticated approach posits the possibility that photodissociation (by 
the ejection of an acetyleneic group) is countered by repair through accretion of carbon 
atoms (Allain, Leach & Sedlmayr 1996a,b). In this case, if Tdi SS is the radiative dissociation 
timescale, and r acc is the C atom accretion timescale then 

Tdiss = (^FUVCdiss) 1 , (15) 

and 

r acc = (nuXck^y 1 , (16) 

where -F FU v is the far-UV radiation field, cr^ss is the effective photodissociation cross section 
per PAH molecule for this particular radiation field spectrum, is the number density of 
hydrogen atoms, Xq is the abundance of C in the ISM, and k acc is the reaction rate for 
sticking of a carbon atom onto such a molecule. For the far-UV radiation field, we take the 
total radiation field above the ionization potential adopted for these molecules (> 6 eV). 
This provides the condition for PAH destruction (Dopita et al. 2002): 

acc j 

or 

n= F E} jy > Xck^ L ^ (i7) 

cn H cr diss 
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where Ti is the Habing Photodissociation Parameter, defined by analogy with the dimension- 
less Ionization Parameter U used in H II region theory. From observation, values of Ti ~ 0.05 
are observed in regions which still contain PAHs (Allain, Leach & Sedlmayr 1996a). However, 
this is a line of sight average, which provides only an upper limit. A threshold of Ti ~ 0.005 
for the destruction of PAHs ensures that PAHs are excluded from the H II regions, but are 
present in the PDRs of even the compact H II regions. However, similar results are obtained 
by assuming that PAHs are destroyed on a short timescale whenever the spectrum-averaged 
photodissociating field exceeds about ten times the local interstellar field. Both of these 
formulations are available within the MAPPINGS Illq code. In this paper we have opted to 
use the Habing Photodissociation Parameter formulation with Ti ~ 0.005, which effectively 
excludes PAH from surviving within the ionized gas. 

When, locally, the threshold of Ti ~ 0.005 is exceeded, we have assumed that the PAHs 
are destroyed in our models, and in this case we return the Carbon that they contain to 
the gas phase. In starburst galaxies, this condition is sufficient to ensure that the PAHs are 
destroyed in both the ionized and in the warm diffuse media. Since the PAH absorption 
cross section is like that of carbonaceous grains, this means that most of the carriers of the 
2200A absorption are removed from lines of sight which permit the UV stellar continuum 
to escape. As we show below, this provides a natural explanation of why the 2200A feature 
is weak or absent in the UV spectra of AGN or starburst galaxies (Calzetti 2001), provided 
that small organic grains or other carbonaceous types are not present in high abundance. 



5. Grain Size Distribution & Composition 

The grain size distribution and the grain composition are the principal factors which 
determine the wavelength-dependence of the absorption and scattering processes. The grain 
size distribution in the ISM results from the balance between the grain formation and de- 
struction processes. It has usually by a power law over a wide range of sizes, a, (Mathis, 
Rumpl, & Nordsieck 1977); 

dN(a)/da = ka~ a a min < a < a max - (18) 

where a, a m \ n and a max are derived by an empirical fit to the scattering and extinction 
in the local ISM. For the MRN distribution the minimum size is set to a min = 0.005yum for 
carbonaceous grains, and to a min = 0.01/xm for silicaceous grains. For both grain types, 
Omax = 0.25/im. 

A more complex distribution has been suggested on physical grounds (Weingartner & 
Draine 2001a), and this takes into account both silicate and carbon-containing grains with 
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different grain-size distributions. For the carbonaceous grains, the PAHs are treated as 
providing an additional component on a continuous distribution of grain sizes. Even these 
distributions can be approximated by power-laws over a wide range of radii. 

Grain shattering has been shown to lead naturally to the formation of a power-law 
size distribution of grains with a ~ 3.3 (Jones, Tielens & Hollenbach 1996), observationally 
indistinguishable from the MRN value, a ~ 3.5. Such a grain size distribution can also hold 
only between certain limits in size, the smaller size being determined by natural destruction 
processes such as photodestruction, and the upper limit being determined by limits on the 
growth by condensation and sticking. To capture these elements of the physics, we have 
adopted a modified grain shattering profile with the form; 

g-(a/a min )~ 3 

dN(a)/da = ka a j-. ^ . (19) 

This creates a smooth exponential cut-off in terms of the grain mass at both the mini- 
mum and maximum grain size of the distribution. This smooth cut removes any edge effects 
in either the emission or extinction by dust that arise due to the sharp cut-offs in the other 
distributions. Our distribution does not have the extra bump to accommodate the PAHs, 
since they are treated as a separate molecular component, as described above 

The grain size limits for both graphitic and silicaceous grains are set to a min = 0.004/im, 
and a max = 0.25/zm. The somewhat smaller size cutoff than the MRN distribution provides a 
stronger dust re-emission continuum below 25 /xm, which is needed to fit the observed IRAS 
colors of starbursts (Rush, Malkan & Spinoglio 1993). The constant k is determined by the 
total dust-to-gas mass ratio, which is determined by the depletion of the heavy elements 
onto dust. Our models use a solar abundance set (as modified by the latest abundance 
determinations Allende Prieto et al. (2001, 2002); Asplund (2000); Asplund et al. (2000)). 
The depletion factors are those used by Dopita et al. (2000) for starburst and active galaxy 
photoionization modeling and are similar to those found by Jenkins, Jura & Loewenstein 
(1987) and Savage & Sembach (1996) in the local ISM using the UV absorption lines to 
probe various local lines of sight. 

Within the MAPPINGS Illq code, the dust grain size distribution is divided into 80 
bins spaced logarithmically between 0.001-10 /zm. The number of grains of each type in each 
bin is then determined by equation 19. The absorption, scattering and photoelectric heating 
is calculated for each size bin. This is then used to establish the temperature probability 
distribution for computation of the FIR re-emission spectrum with quantum fluctuations. 

When present, the PAHs are set at an abundance which uses llOppm of Carbon. Conse- 
quently, the carbonaceous grain component is present at only rather low abundance; equiv- 
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alent to 69ppm of Carbon. 

5.1. Grain Quantum Fluctuations 

In order to properly model the FIR emission bump, it is very important to take stochastic 
quantum heating of the dust grains into account. In particular, the temperature fluctuations 
of grains serves to provide a population of small grains which are hotter than their equilibrium 
temperature and thus an excess in their IR emission at shorter wavelengths. 

This idea of the temperature fluctuations in small grains can be extended to the very 
small grains or molecules, such as Polycyclic Aromatic Hydrocarbons (PAHs). These are 
much more likely to undergo such a fluctuation, although the possibility of destruction by 
a singe photon is also increased. The point at which a dust grain, which cools through a 
continuous radiative decay, should be considered as a dust molecule, which 'cools' (enters 
lower vibrational states) through discrete vibrational and rotational quantum states is rather 
ill-defined. Here we have taken the simple approach and treated all dust as grains having well- 
defined bulk properties, with the molecular treatment being reserved for the PAH population. 

Temperature fluctuations and stochastic heating processes have been considered by sev- 
eral authors since Greenberg first suggested these ideas (e.g. Duley 1973; Greenberg & Hong 
1974; Purcell 1976; Dwek 1986). The most detailed work on dust and temperature fluc- 
tuations has been a series of papers by Draine & co-authors (Draine & Anderson 1985; 
Guhathakurta & Draine 1989) culminating in a set of papers with Li (Draine & Li 2001; Li 
k Draine 2001). 

The FIR emission treatment used in the MAPPINGS Illq code, whose results are pre- 
sented here is based on the algorithms of Guhathakurta & Draine (1989) (GD89) and Draine 
& Li (2001) (DL01). The code provides an optimized solution of their algorithms to deter- 
mine the dust grain temperature distributions for each grain size according to the both 
the strength and detailed spectrum of the local radiation field. The code then integrates 
the resultant FIR emission of the ensemble of dust grains to provide the local re-emission 
spectrum, which is itself integrated through the model in the outward-only approximation. 

6. The Attenuation by Dust; IR to UV 

If our dust model is really representative of the dust found in starburst galaxies, then 
it should reproduce the wavelength-dependent attenuation law seen in starburst galaxies 
(Calzetti 2001). Here we use the word 'attenuation' to distinguish it from extinction, since 
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Table 2. The Solar Abundance Set (Z Q ) and logarithmic depletion factors log(D) adopted 

for each element. 



Element 


log(Z ) 


log(D) 


H 


0.00 


0.00 


He 


-1.01 


0.00 


C 


-3.59 


-0.52 


N 


-4.20 


-0.22 





-3.34 


-0.22 


Ne 


-3.91 


0.00 


Na 


-5.75 


-0.60 


Mg 


-4.42 


-0.70 


Al 


-5.61 


-1.70 


Si 


-4.49 


-1.00 


s 


-4.79 


-0.22 


CI 


-6.40 


-0.30 


Ar 


-5.44 


0.00 


Ca 


-5.64 


-2.52 


Fe 


-4.55 


-2.00 


Ni 


-5.68 


-1.40 
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it properly describes the wavelength-dependent escape fraction of starlight from the galaxy 
averaged across all possible sight lines rather than the wavelength-dependence of the effective 
extinction optical depth along a single sight line. 

Models of the attenuation show that it can be well-represented by the effect of a fore- 
ground screen (Meurer et al. 1997; Meurer, Heckman & Calzetti 1999). In an earlier paper, 
Fischera et al. (2003) showed that many features of the Calzetti attenuation law for star- 
bursts (Calzetti 2001) could be reproduced by turbulent screen having a log-normal column 
density distribution. The grain model used in that paper was that of Weingartner & Draine 
(2001a), which, being strictly applicable to the local ISM, gave a strong 2200A absorption 
feature. However this feature is not evident in the attenuation curves of either starburst 
galaxies or AGN. 

Here we make the hypothesis that the major carrier of the 2200A absorption feature are 
the PAHs, which are, as we have seen, effectively destroyed within the strong UV radiation 
field environments of the diffuse neutral or ionized phases of the ISM. They do, however, 
persist in the PDRs around the boundaries of dense molecular clouds where they are strongly 
excited to fluoresce in the 'UIB' PAH features in the infra-red. These molecular clouds 
are completely opaque to the passage of optical and UV photons, and so provide only a 
geometrical covering factor to the attenuation at these wavelengths. 

This model offers the potential both to explain the absence of the 2200A absorption 
feature and the great intensity of the UIBs in starburst galaxies. Does it work? 

In Figure 3, we show the wavelength dependence of the optical properties of our dust 
model compared with that of Weingartner & Draine (2001a). It is clearly apparent that the 
2200A absorption feature is strongly suppressed by our grain model, but otherwise that the 
extinction and scattering properties of our grain mix is quite similar to the Weingartner & 
Draine (2001a) results once the offset due to the abundance set is accounted for. 

For comparison, we show the optical data for our grain model for the PDRs. Here the 
2200A feature is very strong, and the IR opacity low. A linear combination of the model 

Fig. 3. — The extinction, absorption, scattering and g-factor (thick solid) of our dust model 
for the diffuse phase of the ISM in which the UV field is suflicienty strong for PAHs to 
be destroyed. This component is assumed to provide the foreground screen attenuation for 
the starburst. The absolute to relative extinction is R v = 3.38. The dotted curves are the 
corresponding values for the local interstellar medium of our galaxy as defined by Weingartner 
& Draine (2001a). The generally lower opacity of our grain model mostly reflects the lower 
abundances associated with the revised solar abundance set. 
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curves in Figure 3 and 4 would produce a outcome which is almost indistinguishable from 
those of Weingartner & Draine (2001a). 

We have used the optical parameters of Figure 3 to construct an extinction curve for 
this grain model. This is shown as the thin solid line on Figure 5. We then used the same 
procedure as described in Fischera et al. (2003) to derive an attenuation curve appropriate 
to a turbulent foreground screen having a log-normal distribution in column densities and 
having an i?v-value = 4.0 consistent with that obtained for the Calzetti curve (Ry = 4.04). 
Our theoretical attenuation curve is shown on Figure 5 as the thick solid line, for comparison 
with the Calzetti empirical attenuation curve for starbursts (dashed line). 

It is clear that our grain model is extremely successful in reproducing the Calzetti 
empirical law to an accuracy which lies within the observational errors in the empirical fit. 
This success relies upon both the destruction of the PAHs in the diffuse phases of the ISM, 
and upon the low abundance of other carbonaceous grain types which was forced upon us 
by the need to satisfy C abundance constraints while at the same time ensuring that the 
strength of the PAH emission bands matched observational constraints. 

For the convenience of observers, we have given this theoretical attenuation law in 
tabular form in Table 3. Data values are given in terms of the variable 

, _ A X -A Y _ ff(A-V) 

eA -i^-£(Erv) (20) 

and the absolute attenuation can be calculated from these values by: 

Ax = (Zx/Rv + 1) Ay (21) 

with R y = 4.0. 



7. Modelling the Starburst 



We assume that the starbust region occupies a roughly spherical region filled with H II 
regions of all possible ages - the exact shape of the starburst region does not alter the results 



Fig. 4. — The extinction, absorption, scattering and g-factor for our dust model where 
HOppm of the carbon is assumed to be condensed in PAH. These curves correspond to our 
assumed grain and PAH population in the PDRs and molecular clouds. Note the extreme 
strength of the 2200A absorption, and the deficit in the absorption and extinction at IR 
wavelengths, relative to the Weingartner & Draine (2001a) grain model (dashed curves). For 
this grain model the absolute to relative extinction is R v = 1.92. 
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Fig. 5. — Attenuation curve for a turbulent dusty foreground screen (thick solid line) derived 
using our dust model. The i?v-value is chosen to be 4.0 consistent with the value of the 
Calzetti curve (Ry = 4.04). The Calzetti attenuation curve is shown as dashed line and is 
identified by name. The initial extinction curve for our dust model is shown as thin solid 
line. 



Table 3: Relative Attenuation E(X - V)/E(B - V) given for R v = 4.0 (A v = 1). 



A[/zm] 


Waveband 


This Work 


Calzetti 


Weingartner & 




or Feature 




(2001) 


Draine (2001) 


0.0912 


Lyoo 


9.164 


11.033 


8.591 


0.1053 




8.542 


9.198 


7.308 


0.1216 


Lya 


7.397 


7.765 


5.885 


0.1550 


C IV 


5.999 


6.062 


4.310 


0.1906 


cm] 


4.954 


5.037 


4.471 


0.2175 


2175 


4.839 


4.415 


5.619 


0.2480 




3.904 


3.790 


4.025 


0.3000 




2.759 


2.862 


2.729 


0.3650 


U 


1.839 


1.896 


1.847 


0.440 


B 


1.000 


1.000 


1.000 


0.458 


V 


0.000 


0.000 


0.000 


0.720 


R 


-1.177 


-1.111 


-1.103 


1.030 


I 


-2.432 


-2.267 


-2.211 


1.239 


J 


-2.905 


-2.720 


-2.626 


1.649 


H 


-3.401 


-3.275 


-3.105 


2.192 


K 


-3.682 


-3.690 


-3.444 
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of the modelling. At any age, these H II regions have a radius determined by both their age 
and the pressure in the ISM, according to Figures 1 and 2. The region around the starburst 
is assumed to be filled with a warm, dusty PAH-free diffuse ISM which provides a foreground 
screen extinction as described in the previous section. The starburst is assumed to form stars 
at a rate of 1 M Q yr _1 for a total period of 10 s years, and since each cluster has 10 4 M Q of 
young stars, there are 1000 H II regions younger than 10 Myr. After 10 Myr, the ionizing 
photon output of the central cluster has dropped to negligible levels. Stars older than lOMyr 
are distributed within the H II regions in proportion to the volumes of the individual H II 
regions. The H II regions are assumed not to be 'leaky' to the EUV photons. That is to say, 
they fill their Stromgren spheres and so absorb all the ionizing photons. 

The input stellar spectra (stellar cluster + old stars) within the volume of an H II region 
of given age are computed using the STARBURST99 code as described above. The MAP- 
PINGS Illq code is used with these stellar spectra to compute the radiative transfer through 
the ionized gas, the equilibrium ionization and temperature, and the dust and nebular line 
and continuum spectra throughout the individual H II regions. The MAPPINGS Illq models 
have an isobaric density structure, where the internal pressure is defined by our evolutionary 
modelling (see figure 2), and their physical size matches the computed diameter (see figure 
1). The integrated output spectra are weighted by the number of H II regions in each time 
step, to give the flux equivalent to a galaxian mean star-formation rate of 1 M yr _1 . The 
time sampling of the individual models is first at 0.4Myr, then at lMyr and then at lMyr 
age intervals up to lOMyr. 

We should note here that the MAPPINGS Illq code treats the radiative transfer in the 
outward-only approximation, scattered photons being divided into an outward stream and 
a back-scattered stream which is assumed to return across the ionized nebula. 

Following Silva et al. (1998) and Popescu et al. (2000), we consider each H II region 
to be surrounded by an incomplete opaque shell of molecular gas, in pressure equilibrium 
with the ionized plasma, which blocks a fraction of the incident stellar light according to the 
geometrical covering factor. The PDR at the inner boundary of the molecular shell absorbs 
which almost all the incident FUV, UV and optical photons. To ensure that this is the case, 
the models were continued into the PDR until an H I column of log(nni) = 21.5 was reached. 
This corresponds to Ay ~ 3mag - equivalent to nearly lOmag at the Lyman Limit. Thus, 
effectively all the stellar photons and photons of nebular origin which can excite the PAHs, 
or produce FIR continuum shortward of 100/im are fully absorbed in the PDR. 

If the incomplete shell covers a fraction / = Q/4n of the total solid angle, then in the 
remaining 1 — / of the solid angle, both the remaining (unabsorbed) stellar photons, and 
the photons produced in the H II region are assumed to escape from the starburst region. 
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Some of these will be lost in the outer foreground turbulent dusty screen, but the additional 
attenuation produced here can be calculated using the results of Section 6. The factor 1 — /is 
therefore the escape probability for the stellar photons which successfully exit the starburst 
region. In this sense, the molecular clouds serve to provide a 'grey' attenuation factor of 
1-/. 

Since the hot stars are born in molecular clouds, and these clouds are dissipated during 
the lifetime of the H II region as the region evolves from compact to classical shell morphology, 
it is clear that the molecular cloud covering factor / is a decreasing function of time. This 
concept of age-selective obscuration was first introduced by Silva et al. (1998), and employed 
by Tuffs et al. (2004) in their calculation of the attenuation of starlight in normal galaxies. 
Here, we have investigated two cases. In the first, the covering factor decreases linearly with 
time from 0.9 at 1 Myr to 0.1 at 9 Myr, is complete at zero Myr and equal to 0.1 at times 
later than 9 Myr. Since this function is a little arbitrary, we also investigated the more 
general case in which the covering factor decreases exponentially with time: 

/(f) = exp [-t/r cleax ] (22) 

with the molecular cloud clearing timescale, r c i ea r, set to 1, 2, 4, 8, 16 and 32 Myr. It turns 
out that the linear case provides an SED which is indistinguishable from the exponential 
case and has an effective timescale r ~ 7Myr. 

Since, in general, there will be no physically a priori way of constraining the molecular 
cloud clearing timescale, we are forced to leave it as a free parameter, which can be hopefully 
constrained by the observed SED or colors of real starburst regions. 

To show the way in which the final SEDs are composed of the sums of many H II region 
models, and to demonstrate that all models provide important contributions to the final 
SED, we display the cumulative SEDs resulting from this addition for H II regions of ages 
0.4 - 10 Myr in figures 6 and 7. These figures apply to the case where the molecular cloud 
covering factor decreases linearly with time, and for an ISM pressure P/k = 10 4 cm _3 K. The 
y-axis is chosen to display the energy flux uF u , rather than F u , in order to better show how 
the FIR bump re-distributes the UV stellar continuum. The absolute scale for the summed 
SED corresponds to the star formation rate of 1.0 M yr _1 used in the models. 

Fig. 6. — The cumulative SED for the H II regions and their stars of different ages 0.4 (lower 
curves)- 10 Myr (top). The older stars are included within the volume of the H II regions. 
Each contributing spectrum is weighted in proportion to its covering factor 1 — f(t) The case 
of the covering factor decreasing linearly with time and the ISM pressure P/k = 10 4 cm~ 3 K 
is shown here (see text). 
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In figure 6, which shows the fraction of flux escaping from the H II regions, note the 
strength of the nebular lines and continuum, the weakness or absence of PAH features, the 
hot dust temperatures and the presence of silicate absorption in the infrared spectra of the 
youngest H II region models. For these models, the stellar continuum dominates only in the 
UV. In the visible and near-IR the nebular continuum dominates, with the bound-free edges 
prominent. 

After about 3Myr, the red supergiant stars appear with their characteristic CO features 
in the H and K bands. For these models of older clusters, the stellar continuum dominates 
the SED up to about 5/zm . 

In figure 7, the additional contribution to the SED by the H II region in the remaining 
solid angle and by the photodissociation regions is shown. In these spectra, all the stellar the 
UV continuum radiation has been absorbed in the PDR and converted to dust re-emission, 
which thanks to the low optical depth in the IR is free to escape. Note how the peak of 
the FIR dust re-emission bump shifts to longer wavelengths for the models with the older 
central clusters. This is a consequence of both the decrease in the stellar luminosity and 
the larger shell sizes, which both lead to lower dust temperatures. Also, the PAH emission 
bands are stronger in the older H II regions. In these models, the PAHs survive thanks to 
the lower photo-destruction rates of the PAHs resulting from the lower UV fluxes and cooler 
stellar temperatures. Nonetheless, there are still plenty of photons with 6 < hv < 13.6eV to 
be absorbed by the PAHs and to excite the UIB features. 

Thus our models can succeed in explaining the apparently contradictory observations 
that starburst galaxies are seen to have very strong UIB emission (implying the presence 
of large quantities of photon-excited PAH molecules), but have no sign of the UV 2200A 
absorption feature (implying the absence of PAHs, or at least carbonaceous grains, in the 
more diffuse phases of the ISM). 

Fig. 7. — The cumulative SED of figure 6 (lower curve) to which has been added the contri- 
butions of the optically thick portions of the H II regions and their PDRS with ages from 0.4 
(lower curves) to 10 Myr (top). As in figure 6, the case of the ISM pressure P/k = 10 4 cm _3 K 
is shown, and the case for which the covering factor changes linearly with time (see text). 
Each contributing spectrum is added with a weight in proportion to its covering factor fit). 



-29- 



8. Parameters Controlling the SED 

In this section, we present the results of the systematic investigation of the two param- 
eters which, apart from the optical depth of the foreground absorbing screen, the assumed 
dust grain size distribution and the assumed dust composition, primarily control the SED 
of the starburst. These parameters are identified as the mean pressure in the ISM, and the 
molecular cloud covering factor. 

8.1. The Role of Pressure 

We have computed three sets of models with different pressures, P/k = 10 4 , 10 6 and 
10 7 cm~ 3 K. The first of these corresponds roughly to the pressure in the ISM near our sun 
(Jenkins 1983), and so would be applicable to star formation complexes in disk galaxies. 

Pressures of P/k ~ 10 6 are commonly encountered in starburst galaxies, and can be 
measured using either the emission measure in the H II regions or even using the density- 
sensitive [S II] AA6717/673lA ratio, provided that n e > 100 cm~ 3 . Starbursts with these 
pressures are characterized by warmer IRAS colors (Kewley et al. 2001, b). 

Finally, starbursts with P/k ~ 10 7 have the warmest IRAS colors, and are most often 
the nuclear starbursts in post-merger galaxies. For these, the [S II] AA6717/673lA ratio 
indicates n e > 1000 cm -3 (Kewley et al. 2001, b). Thus, our choice of parameters should 
cover the full range seen in real starburst galaxies. 

Figure 8 shows the effect of the changing pressure on the SEDs. The main change is the 
progressive shift of the peak of the FIR dust re-emission bump to shorter wavelengths, and 
the increase in its flux. For a pressure of P/k = 10 4 , the peak is at 105/im, at P/k = 10 5 , 
the peak is at 70/im and by P/k = 10 7 , the peak has shifted down to 50fim. At this high 
pressure, hot dust emission can be traced to as wavelengths as short as 4/xm. 

For all these curves, there is no discernable change in the stellar continuum. This is to 
be expected as it is defined by the geometrical constraints. However, there is a progressive 
change in the nebular line spectrum towards higher excitation at higher pressure. This 
is a consequence of the mean ionization parameter, which slowly increases as the pressure 

Fig. 8. — The summed starburst SEDs computed for ISM pressures P/k = 10 4 , 10 6 and 
10 7 cm~ 3 K. Computed for the case in which the molecular cloud covering factor changes 
linearly with time. 
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increases. 

Surprisingly enough, the strength of the PAH emission features is almost invariant with 
pressure. This arises because the PAH molecules can survive in the PDRs independently 
of the pressure, since the photodissociation rates are low in these zones. Thus, the PAH 
features scale only with the covering factor, /. 

We would expect these theoretical SEDs to be not very accurate at wavelengths > 100/im 
because we have not taken into account either the dust re-emission of diffuse starlight, or 
caused by stars older than 10 8 years. Nor have we included the FIR re-emission by foreground 
dust outside the main star-forming region. Thus, we would expect that in real starburst 
galaxies there will be a cold component of dust emission from the ambient medium which is 
not predicted by this model. However, for starburst galaxies the results in Section 9 shows 
the fit to be remarkably good. Nonetheless, the absence of this component will certainly 
become a serious problem for the modelling of "normal" (non-starburst) galaxies, since we 
know that for these galaxies the luminosity of the cold dust component generally exceeds 
that of the locally heated dust e.g. Popescu et al. (2002). High resolution maps of galaxies 
like M31 (Haas et al. 1998), M33 (Hippelein et al. 2003) or NGC891 (Popescu et al. 2004) 
made in the ISO 170 or 200/xm bands all show that most of this cold dust luminosity is 
carried by the diffuse ISM, not by the optically thick molecular clouds. Thus, the model 
presented in this paper is not really applicable to normal galaxies. For this, the output SED 
of the HII regions computed here will have to be incorporated in a model such as that of 
Popescu et al. (2000) that self-consistently calculates the radiation transfer in the diffuse 
component of the ISM. 

Another way in which these theoretical SEDs may differ from those of real starburst 
galaxies is in the the way collective effects may alter the form of the computed SED. In 
galaxies with a strong centrally concentrated starburst, H II regions may be superposed 
along the line of sight, and the radiative transfer would be then more complex than assumed 
in our model. Again, the solution to this is to build a three-dimensional radiative transfer 
model for the dusty gas, and to use the SEDs presented here as input source functions. 

Finally, there may well be a population of dust-enshrouded and confined ultra-compact 
H II regions contributing to hot dust emission. We might reasonably expect to find relatively 
more of these in high-pressure regions. In our own Galaxy the W3 complex provides pre- 
cisely these conditions. At high pressure, the star-forming regions are embedded in denser 
molecular gas dust, and the H II regions stall sooner, and at smaller radii. Single OB stars 
would therefore find it difficult to blow their H II regions to sufficient size to be able to 
communicate with regions of lower density, and may therefore spend a significant part of 
their main-sequence lifetime in the Compact or Ultra-Compact H II region phases. We will 
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return to the evidence for such a population of compact H II regions in Section 9, below. 

8.2. The Role of Covering Factor 

The molecular cloud covering factor has a significant effect on nearly all the starburst 
SED, as Figures 9, 10, and 11 make clear. These figures show the effect of the exponential 
molecular gas clearing timescale, r clcar at the three different pressures. 

First, high covering factors lead to an attenuation of all the stellar continuum, which 
is dominant from the Lyman Limit down to about 5/im. However, when the molecular gas 
curtain is slow to open, there is a significant effect of ageing in the stellar continuum, since 
the light from the oldest H II regions is dominant, and these H II regions, being largest, also 
contain the largest number of stars with ages > lOMyr. This ageing is visible as an increase 
in the size of the Balmer and Paschen discontinuities, and a flattening of the UV spectral 
slope. 

Larger covering factors also lead to an increase in strength of the PAH features, and both 
a strengthening and broadening of the ~ 100/im dust emission peak. This strengthening of 
the peak works in the opposite direction from pressure, since the strengthening is due to an 
increasing contribution of the older H II regions with longer r c i car . These older H II regions 
are larger, and consequently the dust in them is cooler. 

At about 25/im, the flux changes little with r c i ea r for all pressures. However, the flux at 
this wavelength is very strongly affected by affected by changes in pressures. By contrast, 
the stellar continuum and the PAH features are strongly affected by the molecular cloud 
covering factor, but changes in pressure do not affect these at all. This offers the potential 
to distinguish observationally between these two parameters, as we will discuss below. 

8.3. Approximations to the Nonthermal Emission 

One of the results of our modelling is that the 60/im flux depends quite strongly on 
the pressure in the interstellar medium. However, we also know from observation that there 

Fig. 9. — The summed starburst SEDs computed for ISM pressure P/k = 10 4 cm~ 3 K with 
exponential molecular cloud clearing timescales of 1, 2, 4, 8, 16 and 32 Myr. Note that a 
long clearing timescale depresses the UV continuum but leads to stronger PAH emission as 
well as stronger ~ 100/im dust emission. 
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is an extraordinarily close correlation between the 60/zm infrared continuum and the radio 
1.4 GHz continuum of star forming galaxies. This linear correlation spans ~ 5 decades 
of magnitude with less than 0.3 dex dispersion (Yun, Reddy & Condon 2001; Wunderlich, 
Wielebinski, & Klein 1987), and the IR -radio correlation persists even when measurements 
of the cold dust emission (which carries the bulk of the FIR luminosity for normal galaxies) 
are included in the total FIR luminosity (Pierini et al. 2003). The mean relationship between 
the 60/zm flux and the 1.4 GHz continuum is: 

LlAGRz 

W Hz" 1 

We know that, in star forming galaxies at 1.4 GHz, the non-thermal emission by rel- 
ativists electrons dominates by at least an order of magnitude over the free-free emission 
(Condon 1992; Niklas, Klein & Wielebinski 1997; Dopita et al. 2002). Therefore, and some- 
what remarkably, this relationship couples a purely thermal process with a non-thermal 
process, over many decades of flux, and locally, within individual galaxies. Because of the 
intrinsic interest of this radio : FIR correlation, and because we have shown that ISM pres- 
sure directly affects one of the parameters that enters into this correlation, in this section 
we will investigate a 'toy' model of the synchrotron emission in order to gain insight as to 
how the synchrotron emission might be similarly affected by the ISM pressure, thus keeping 
the radio : FIR correlation as tight as is observed. 

If the lifetime of the synchrotron electrons is short compared with the evolution timescale 
of the starburst, then as Bressan, Silva & Granato (2002) have shown, the synchrotron 
emissivity acts as a bolometer of the supernova rate, since the relativistic luminosity is 
ultimately derived from the Fermi acceleration process in supernova shocks. Specifically, the 
non-thermal luminosity is given by 

L NT (v) oc M,5 Q -V- a (24) 

where M* is the total star formation rate. Since the observed spectral index lies in the range, 
0.5 > a > 1.0, there is only a very weak dependence on the magnetic field in this case. 

If the lifetime of the relativistic electrons is long compared with the length of the star- 
burst, then the emissivity of the relativistic electrons depends on the square of the local 
magnetic field, and the local density of the relativistic electrons. Groves et al. (2003) as- 
cribes the overall radio-FIR correlation, and its local variation in individual galaxies to the 
operation of magneto-hydrodynamic turbulence. 

Fig. 10. — As Figure 9, but with P/k = 10 6 cm" 3 K. Note the displacement of the FIR bump 
towards shorter wavelengths. 
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In a magnetized and turbulent interstellar medium, numerical simulations (e.g. Cho & 
Vishniac (2000b)) have shown that the local magnetic field will stay in equipartition with 
the gas turbulent energy; 

8V ~ B/^/Arp. (25) 

This provides the coupling between the gas density and the magnetic field which is a 
necessary, but not a sufficient, condition for the operation of the radio - FIR correlation. This 
equation also ensures that the magnetic field pressure and the gas pressure are correlated. 

The local synchrotron emissivity at frequency v is given from the standard synchrotron 
theory by 

j v = f{a)kB ia+1 ^ 2 v- {a - 1 ^ 2 (26) 

where f(a) is a complex numerical factor of a, the power-law slope of the non-thermal electron 
the number density of the relativistic electrons with relativistic 7; N(^)' a = k^~ a , and B 
is the local magnetic field. The total density of the relativistic electrons is therefore n cr = 
^7mTn/ ( a — 1); where 2 < a < 3 corresponds to a frequency spectral index of synchrotron 
emission 0.5 > a > 1.0. Provided that the relativistic electrons cannot escape from the star- 
forming region, the constant, k, should be proportional to the star formation rate. Equation 
26 predicts, by contrast with equation 24, there should be a strong correlation between the 
emissivity and the magnetic field, in the range B 3 ^ 2 to B 2 . 

At normal ISM pressures, the lifetime of the relativistic electrons should be longer than 
the star formation timescale, so that equation 26 would apply. However, as the pressure 
increases, the emissivity of the non-thermal emission also increases, and the lifetime of the 
electrons becomes shorter. Thus, at high P/k, we expect that equation 26 ceases to apply, 
and we go over to equation 24 instead. Unfortunately we have no way of estimating at which 
point this occurs with this simple theory. 

In our model for the emissivity, we have assumed that equation 26 applies, and have 
adopted the case a = 2, which provides a spectral index of a = 0.5 and a pressure dependence 
of the emissivity arising from the B— field of P 3 / 4 . The calibration of the non-thermal flux 
was accomplished empirically by applying equation 23 at a pressure of P/k = 10 4 cm~ 3 K. 
The emissivity was then computed for the higher pressures, still assuming that equation 26 
applies. 

The result is shown in figure 12. Here we have plotted the computed F u from the Lyman 
Limit up to 30cm (0.98GHz). The apparent change in spectral index of the radio continuum 

Fig. 11.— As Figure 9, but with P/k = 10 7 cm~ 3 K. 
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with pressure is due to the underlying contribution of the radio free-free continuum from the 
ionized plasma in the H II regions. This is relatively more important at low ISM pressure. 
Changes in the radio non-thermal luminosity with pressure are correlated with changes in the 
60/im flux. However, it is clear that the size of the increase in 1.4GHz continuum computed 
using equation 26 is appreciably larger than the corresponding increase in the 60/im flux. It 
seems clear that when the pressure of the ISM reaches P/k ~ 10 6 cm _3 K, the lifetime of the 
synchrotron electrons is becoming, or has become, shorter than the age of the starburst. 

Once this condition is reached, then the synchrotron flux becomes more sensitive to 
the instantaneous injection rate of relativistic electrons. We would therefore expect more 
dispersion in the radio - farlR correlation at the highest pressures. Since many of the ULIRG 
galaxies have high pressure, this could contribute to scatter at the upper end of the radio - 
farlR correlation c.f. Bell (2003). 

9. Comparison With Observations 

9.1. Pan-spectral fits to Starburst Galaxies 

In this section we compare the theoretical SEDs with the observations of two exceedingly 
well-observed objects; Arp 220 (z = 0.0182) and NGC 6240 (z = 0.0245). These objects are 
chosen because the SEDs of these objects have been established from the Lyman Limit up to 
beyond 1000/im. Arp 220 is widely used as the template for violent and highly extinguished 
starbursts and for the computation of the fluxes expected from high redshift objects. The 
case of NGC 6240 (IRAS 16504+0228) is not so clear-cut. It is is a very well-known ultra- 
luminous infrared merger remnant with a strong non-thermal radio excess. It is clearly 
dominated by a starburst, and shows a very strong Balmer lines but has a LINER or Sy2 
nuclear spectrum (Kewley et al. 2001b; Goldader et al. 1997; Kim et al. 1995; Veilleux et al. 
1995). It also shows evidence for an obscured AGN at X-ray and radio wavelengths (Ikebe 
et al. 2000; Risaliti et al. 2000; Kewley et al. 2000). 

For both galaxies the UV and optical fluxes were taken from the third reference catalogue 

Fig. 12. — The summed SEDs (F v vs. A) computed for ISM pressures P/k = 10 4 , 10 6 and 
10 7 cm _3 K for the case in which the molecular cloud covering factor changes linearly with 
time and for an absolute star formation rate of 1.0M o yr _1 . Note that changes in the radio 
non-thermal luminosity with pressure are positively correlated with changes in the 60/im 
flux. 
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of bright galaxies v3.9 (de Vaucouleurs et al. 1991). The near-IR (JHK-b&nd) points are 
from Spinoglio et al. (1995), and include aperture corrections to allow a direct comparison 
with the larger aperture mid- and far-IR fluxes. The L-band magnitude for NGC 6240 (Allen 
1976) has been similarly corrected by scaling with the same factor as required for the K- 
band. The 3-1500 fim data points were taken from Klaas et al. (1997, 2001), Benford (1999), 
Spoon et al. (2004) and references therein. All UV to near-IR fluxes have been corrected 
for Galactic extinction using the E(B-V) values based upon IRAS 100 /im cirrus emission 
maps (Schlegel, Finkbeiner, & Davis 1998) and extrapolating following Cardelli, Clayton, & 
Mathis (1989). 

The process of matching the models to the observations required two steps. First, the 
data were shifted to rest-frame frequencies and the absolute flux determined from the lumi- 
nosity distance of the galaxy was scaled to the far-IR section of the SED. The location of the 
peak and the slope of the long- wavelength tail of the thermal dust re-emission fixes the ap- 
propriate value of Pj k. The molecular cloud clearing timescale r c i ear was best determined by 
the fit in the 2-5 /im region, which is dominated by the stellar continuum leaking out between 
the molecular clouds. The choice of clearing timescale also changes the inferred scaling factor 
in the far-IR, but to a limited degree. The final scaling factor to be applied to the observed 
absolute fluxes is a direct measure of the star formation rate in the galaxy concerned, since 
all the theoretical spectra are scaled to an absolute star formation rate of l.OM yr _1 . To 
estimate this scaling factor we assume a flat universe with H = 71 km s -1 Mpc -1 , = 0.27 
and Q A = 0.73 (Spergel et al. 2003; Tonry et al. 2003). 

Second, when a best-fit match was obtained throughout the far-IR regime, the model 
SEDs were reddened according to the foreground turbulent screen attenuation law given in 
table 3. The total attenuation was adjusted to match the integrated fluxes at optical and 
near-IR wavelengths. This provides a direct estimate of Ay to the star forming region. 

The family of SEDs generated for this paper and used in this section are available as tab 
delimited text files in the electronic version of the paper. Likewise, the attenuation curves 
from section 6 are also given in the electronic version of the journal in the same form and over 
the full range of wavelength so as to facilitate their use by others. In these files, and in the 
fitting used in this section, we have not included any synchrotron contributions, since these 
are so uncertain at present. Without this contribution, the model SEDs become dominated 
by the free-free emission (assumed optically thin) at wavelengths longer than about 2 mm. 

The observed SEDs of Arp 220 and NGC 6240 with model fits overlaid are shown 
in figures 13 and figures 14. High pressure models with P/k = 10 7 cm~ 3 K were required 
in both cases. The SED of NGC 6240 could be accurately fitted with either a linearly 
decreasing covering factor or r c i car ~ 8 Myr and with Ay = 2.7. However, for Arp 220 a 
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good match to the far-IR could be only be obtained for a model with r c i ear = 100 Myr. This 
implies that the starburst is enshrouded by an optically thick dust and molecular screen 
which absorbs almost all UV/optical radiation. This molecular screen effectively acts as a 
bolometer. The fit in the near-IR regime requires an additional attenuation of Ay = 5.7. 
It was not possible to fit the excess of flux in the UV. This component may arise from an 
outer stellar population which is not as deeply dust embedded as the main star formation 
region. If we use the value of N H /A B — Ay — 5.8 x 10 21 cm~ 2 mag -1 (Bohlin et al. 1978) 
along with Ry = 4.0, then our value of Ay ~ 5.7 mag translates to an inferred column 
density of Ah ~9x 10 21 cm~ 2 . This can be directly compared with the X-ray absorption 
to the nuclear source, N H ~ 3 x 10 22 cm -2 (Clements et al. 2002). Since we would expect 
the central source to be more deeply dust-embedded than the more extended star formation 
region, our estimate of Ay = 5.7. 

With luminosity distances of 78 Mpc and 105 Mpc for Arp 220 and NGC 6240 re- 
spectively, we infer star formation rates of 300 M yr _1 and for Arp 220 and 120 M Q yr _1 for 
NGC 6240. As one might expect for these heavily dust-enshrouded objects, these are very 
close to the estimates derived from the bolo metric far-IR luminosities, Lir. 

The fit with the observations is, in general, rather good. For Arp 220, the main differ- 
ences between the the data and the models are the excess in UV flux in the data, already 
alluded to, and a difference in the observed and predicted strength of the PAH features. The 
data show that the strength of the PAH features are overestimated by a factor of about two 
in the models. This difference can be understood as a result of the very high star formation 
rate inferred for Arp 220, which implies a strong PAH- photodissociating UV field in the 
star formation region. Presumably, a greater fraction of the PAHs are destroyed in this 
environment. 

Relative to the model, Arp 220 also shows an excess of flux in the 12 — 30/im band 
which can be ascribed to the presence of an additional warm dust component. As discussed 
in Section 8.1, this is most likely due to a population of confined compact and ultra-compact 
H II regions around single OB stars, or small clusters. It is highly unlikely that this is due 
to the presence of an AGN, since at other wavelengths, there is no evidence to support this 
idea. 

Fig. 13. — The global fit to the spectrum of Arp 220. This requires a P/k = 10 7 cm" 3 K, 
a very long clearing timescale for the molecular clouds, and a foreground dust screen with 
Ay ~ 5.7 mag. The fit implies an absolute star formation rate of 300 M yr _1 . Also shown 
is the modified Black Body fit to the 60 — 500 /xm region. 
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For NGC 6240, the main difference between model and observation is that the obser- 
vations also show an excess of flux in the 10 — 30/zm waveband. In this object, however, an 
AGN is undoubtedly present. 

Overall, the goodness of fit over the whole range of wavelength was somewhat surprising 
to us. In particular, we would have expected to find a poor quality of fit in the 100 — 1000/im 
region of the spectrum, since here the emission was thought to be dominated by the cooler 
dust located in the more extended dusty region surrounding the starburst, which is not 
represented in our modelling. However, this can be understood as a consequence of the 
long molecular cloud clearing timescales we inferred in the fitting procedure. This implies 
that the vast majority of the UV photons are being absorbed in the PDRs surrounding the 
H II regions. This is consistent with the rather long molecular cloud clearing timescales 
inferred for these objects. Because of the high pressures these regions are very dense, and 
consequently, thin in comparison with the radius of the H II regions. Thus the intensity 
of the local radiation field is simply a function of the optical depth into the PDR. Since 
the local radiation field determines the local temperature, or more strictly, the temperature 
distribution of the dust grains, the global dust emission spectrum is a superposition of a wide 
range of dust temperatures. The hottest grains are located close to, or within, the ionised 
region, and the coolest grains are in the tail of the PDR. 

Frequently, observationalists use a single temperature modified blackbody spectrum as a 
fit to the 100 — 1000/im spectra of starbursts. We have therefore generated such spectra with 
dust emissivity proportional to A _/3 . These have been fitted over the region 60 — 500/xm. 
Both fits yield a value (5 — 1.5 (similar to what is usually assumed). The effective dust 
temperatures derived are T = 37.5 K for Arp 220 and T = 41.3 K for NGC 6240. This is in 
good agreement with what has been found for single temperature fits to bright IRAS galaxies 
(Dunne et al. 2000). Estimates of the Bolometric flux based on such parameterizations 
clearly underestimate the true value, since these fits do not take account of the emission at 
the shorter wavelengths. More complex multi-temperature fits do a better job in this regard 
(see e.g. Klaas et al. 2001), but all such fitting is fundamentally unphysical. Especially 
misleading is the use of an additional parameter, the optical depth at 100 /xm. If this is 
mis-interpreted as representing a real physical parameter, it would imply such tremendous 
optical depths that the starbursts would be totally invisible at optical wavelengths. Such is 

Fig. 14. — The global fit to the spectrum of NGC 6240. This requires a P/k = 10 7 cm _3 K, 
a molecular clearing timescale r c i car ~ 8 Myr and Ay = 2.7. The fit implies an absolute star 
formation rate of 120 M yr _1 . Also shown is the modified Black Body fit to the 60 — 500 /mi 
region. 
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clearly not the case. 

Nonetheless, the single-temperature fits are a useful way of characterizing the form 
of the long wavelength side of the far-IR thermal emission from dust. For example Blain 
et al. (2004) have used a sample of ultra-luminous galaxies in the nearby universe as well 
as a sample of more distant starburst galaxies (1 < z < 5) to find an (admittedly loose) 
correlation between the inferred dust temperatures and the luminosity of the host galaxy. 
In the framework of the theory presented here, this could be interpreted as a correlation 
between the star formation rate and the pressure in the ISM of the host galaxies. 

9.2. IR Color-Color Diagrams for Starburst Galaxies 

In the following sections we compute theoretical color : color diagrams for our model 
starbursts as would be seen using IRAS or the Spitzer Space Observatory. We regard the 
computed color: color diagrams involving flux densities at wavelengths longer than 100/im as 
unreliable such as those involving the Spitzer Space Observatory MIPS 160/im filter). This 
is because we have neglected both the old stellar populations and the diffuse components 
of the dust, both of which are important contributors at these wavelengths for starbursts 
located in more normal galaxies (Popescu et al. 2000). 

9.2.1. IRAS 

We have used the band pass data for the IRAS 12, 25, 60 and 100/xm bands in conjunc- 
tion with our theoretical SEDs to construct color-color diagrams. These are shown in figure 
15 . 

Some homogeneous observational data is shown for comparison. We have used Rush, 
Malkan & Spinoglio (1993) data for objects classified as starbursts. These are plotted as 
symbols. Typical errors in this data set amount to 0.1 to 0.2dex in terms of errors on 
color-color plots. 

The single square represents the mean for the Rush, Malkan & Spinoglio (1993) Seyfert I 
colors as determined by Dopita et al. (1998). It is clear that our theoretical curves define 
an upper envelope to the observed points on all four of these color-color diagrams. However, 
our detailed spectral fitting in Section 9.1 shows that our models are missing a warm-dust 
emission component in the 12 — 30 /zm region, which is almost certainly caused by a popu- 
lation of confined compact and ultra-compact HII regions with ages less than about 2 Myr 
surrounding single OB stars or small OB clusters. This could account for the offset between 
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the theoretical curves and the majority of the Rush, Malkan & Spinoglio (1993) starbursts. 
However, some of the objects classified as starbursts by Rush, Malkan & Spinoglio (1993) 
are in fact of mixed excitation; starburst + AGN. Indeed, the fact that those objects that 
lie particularly low on one of these color-color plots also lie low on the others makes this 
conclusion secure. From the mixing lines given in Dopita et al. (2004) we might infer that 
these objects have a bolometric contribution ascribable to the AGN of between 10 and 20% 
that of the starburst. 

It is clear that the models reproduce the intrinsic spread in colors observed in real 
starburst galaxies. In general, warmer IRAS colors correspond to higher ISM pressures, and 
therefore we might expect to find a correspondence between the compactness of the starburst 
and its IRAS colors. However, none of these diagrams is suitable for separating the two 
controlling parameters of the ISM pressure and the molecular cloud clearing timescale. 

The scatter of the points on the -F100/-F25 vs ^60/^25 diagram is low. However, this can 
now be seen as simply a degeneracy in the parameters - both P and r c i ear spread points in 
the same direction, which also happens to lie in the same direction as the AGN : starburst 
mixing line. The same also applies to the F W0 /F 12 vs F 60 /F 2 $ diagram. The remaining two 
color-color diagrams offer better potential to separate the contribution of any AGN to the 
total spectrum, and also the pressure can be somewhat constrained. 

9.2.2. Spitzer Observatory 

The IRAC and MIPS instruments of the Spitzer Observatory offer the potential of 
observationally separating and measuring the both the pressure and the molecular cloud 
covering factor, and also of estimating the contribution of any AGN component. 

This can be done in the following manner. First, because the IRAC 4.5/im band directly 
measures the stellar continuum, the 5.7/zm band sees the PAH emission and the 7.9/zm band 
is completely dominated by the PAH emission, the ratio of fluxes in these wavebands is very 

Fig. 15. — Color-color diagrams computed for the IRAS bands. The observational points 
come from the Rush, Malkan & Spinoglio (1993) starbursts. The single point shows the mean 
colors of Seyfert I galaxies given by Dopita et al. (1998). Note that our models reproduce 
the intrinsic range of FIR colors seen in starburst galaxies, and that the grid of models 
in all cases defines an upper boundary around the observed points. It is clear that, if the 
theoretical grid is correct for "pure" starbursts, then a number of the objects classified as 
starbursts are in fact of mixed excitation, starburst and AGN. 
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sensitive to the covering factor. 

Second, the MIPS 24fim band is almost invariant with covering factor, but is extremely 
sensitive to pressure (see Figure 9 through 11). Thus, the flux in this waveband can be used 
in conjunction with either any of the IRAC bands to provide a pressure discriminant. In 
practice either the IRAC 4.5yum or the 5.7/im bands are best for this purpose. Alternatively, 
longer-wavelength MIPS bands could be used, but these tend to suffer the same degeneracy 
as the IRAS bands. 

Finally, the presence of an AGN continuum would tend to wash out the PAH features 
(Laurent et al. 2000), so the ratios of fluxes in the IRAC 7.9 /zm and 4.5/im bands and any 
of the IRAC bands could be used to help determine the AGN fraction. 

Some useful Color-color plots are shown in Figure 16. 

This figure shows how either the F5.7/F4.5 vs -F24/F4.5 or the -F7.9/F4.5 vs -F24/F4.5 color- 
color diagrams provide a clean separation between P and r clcar , while the MIPS color-color 
diagrams, F 70 /F 160 vs F^q/F^ or -F 160 /F 24 vs F 70 /Fi 60 are essential degenerate in these 
parameters. However, by providing a starburst line on these diagrams, they can be used to 
help estimate the contribution of an AGN, once the mean AGN colors have been determined. 



10. Calibration of Star Formation Rates 

10.1. Star Formation Rates from Ha 

Neglecting the internal extinction by dust within the H II regions, the Ha flux can be 
fit in terms of the molecular cloud dissipation timescale. From our models, this gives, to 
better than 1% for r > 2 Myr: 



log 



ma 



erg.s" 



41.63 + log 
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M yr- 1 



log[1.6 + 



r 



Myr 



(27) 



Fig. 16. — Color-color diagrams computed for the IRAC and MIPS instruments on the 
Spitzer Observatory. In general the Spitzer wavebands are more useful in separating the 
effects of ISM pressure from the time-dependent molecular cloud covering factor. Increasing 
pressure increases the MIPS 70 and 24 /iin fluxes, but leaves the IRAC bands unchanged, 
while increasing covering factor decreases the stellar continuum (the IRAC 4.5/zm flux) but 
increases the strength of the PAH bands, most clearly in the IRAC 7.9yum band. 
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and for r = 1 Myr: 
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For the case r = 1 Myr this gives the following calibration for the star formation rate, 
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which can be compared directly with the earlier calibrations by Dopita & Ryder (1994); 
Kennicutt (1998); Panuzzo et al. (2003) (c.f. equation 2). The agreement with these earlier 
works is satisfactory, given the different stellar models used here, and difference in other 
assumptions such as the IMF. Here we require a somewhat larger Ha flux to deliver a given 
star formation rate. However, it should be emphasized that for other values of the molecular 
cloud dissipation timescale, equation 27 should be used, and these will require a greater star 
formation rate to produce a given Ha flux. 



10.2. Star Formation Rates from the UV & GALEX 

The GALEX satellite offers two wavebands centered at 1500A (FUV) and 2200A (NUV). 
We have folded the response function of these band passes with our predicted FUV SED. 
As can be seen from figure 8, the SED is unaffected by the ISM pressure. Instead, it is 
strongly influenced by the molecular cloud clearing timescale. This is because, to first order, 
the molecular clouds act to obscure some fraction of the exciting stars. Since, for any 
particular star, this obscuration is either non-existent or total, the attenuation produced 
by the molecular clouds is effectively wavelength independent. The small changes in the 
slope of the FUV spectrum, and in the relative intensity of particular spectral features in 
the UV, apparent in figure 10, are produced by age differences in the observable stellar 
populations. In regions with long molecular cloud clearing timescales, the younger stars 
tend to be preferentially obscured. 

In real starburst galaxies, the FUV spectrum will also be obscured by the overlying 
foreground screen of material. We will assume that this can be described by the attenuation 
curve given in figure 5 and table 3. Here we are concerned only with the conversion of the 
absolute UV flux to star formation rate. 

We find that for the model starbursts, the GALEX 1500A and 2200A bands provide 
nearly the same flux; F(1500)/F(2200) = 1.03 ± 0.015, for all values of the molecular cloud 
dissipation timescale. The F(1500) flux can be fit, like the Ha flux in terms of the molecular 
cloud dissipation timescale. For all computed values of r, this is given to better than 3% 
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accuracy by: 



log [L1500/ erg s x Hz 



28.62 + log [uf^r] - Zo»[2.3 + 0.37 (^)] 



This flux accounts for the internal dust extinction (or, more properly, attenuation) of 
the H II regions, but it does not account for any foreground screen attenuation. 

The dependence of the flux on r is weaker than for Ha because the Ha flux is produced 
mostly in the first 4-5 Myr of a starburst, and so is very prone to obscuration by the placental 
molecular clouds of the starburst itself. By contrast, the stellar UV flux at 1500A remains 
appreciable out to an age of tens of Myr and so is much less affected. Indeed, because of this 
the computed F(1500) flux is somewhat sensitive to the assumed duration of the starburst, 
with older stars contributing significantly to the UV flux at this wavelength. For example, 
with t = lMyr, a 1.0 M Q yr _1 starburst continued for 10 8 yr will give and absolute flux 
F(1500) = 1.58 x 10 28 erg s -1 Hz -1 . However, if the duration of the starburst is only 10 7 yr 
the UV flux falls to F(1500) = 1.04 x 10 28 erg s" 1 Hz" 1 . 

For comparison with earlier calibrations, for a 10 7 yr duration starburst with a molecular 
cloud dissipation timescale, r = 1 Myr: 
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(30) 



Again, we require a larger observed L ° flux to deliver a given star formation rate, 

to ' H to 1500A to ' 

but larger values of the molecular cloud dissipation timescale require a greater star formation 
rate to produce a given L ° flux. 

1 to 1500 A 



10.3. Star Formation Rates from the FIR 



In the far-IR, fluxes are affected by both the pressure and the molecular dissipation 
timescale, so the situation with respect to determining star formation rates is somewhat 
complex. If we were to choose a single wavelength for which the effect of the molecular 
dissipation timescale is largely compensated, then we would choose the IRAS 25/im or the 
Spitzer Observatory MIPS 24/im band. This is because at this wavelength, the reduction in 
the strength of the PAH bands for shorter molecular dissipation timescales is compensated by 
the increasing color temperature of the FIR re-emission continuum (see figures 9, 10 and 11. 
This waveband probably gives the most reliable conversion of flux to star formation rate when 
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r is not known. In table 4 we give the computed luminosities through the MIPS 24/zm band 
and the IRAS 25/xm band for a star formation rate of 1.0 M yr" 1 . With this star formation 
rate these IR luminosities can be fit to an accuracy of ±15%, which is probably sufficient for 
most astronomical purposes, by a simple function of pressure /og[L 2 5 Mm /L Q ] = f(P/k) with 
f(P/k) = 7.30, 7.69 and 8.00 (MIPS) and 7.66, 8.02 and 8.30 (IRAS) for the three pressures 
P/k = 10 4 , 10 6 and 10 s , respectively. 

More usually, we are accustomed to using the conversion between far-IR luminosity I/fir, 
or the total IR luminosity L IR . For the IRAS wavebands these are defined by (Saunders & 
Mirabel 1996) as: 

LfIR = 2.58L60^m + -^lOO^m (31) 

and 

-^ir — 13.48L 12(tim + 5.16L 2 5 (tim + 2.58L 60Mm + -^ioo^m- (32) 

With these definitions, we have computed the model fluxes for a star formation rate of 
1.0 M yr _1 . These are given in table 5. There is no simple fit to these luminosities, owing 
to the sensitivity to r, particularly at the larger values of P/k. 

When r is long, and the pressure is high, then the majority of the stellar flux is re- 
radiated in the far-IR, and Lir becomes a good measure of the bolometric luminosity of 
the exciting stars. This is particularly true at high pressure, since the dust temperature 
is then high enough that the emission beyond 100 /im becomes unimportant. Put another 
way, the dust shell is virtually complete around the stars and so the dust shell acts as an 
IR bolometer. For the most extreme case computed (r = 32Myr; P/k = 10 7 cm~ 3 K), L IR 
provides a conversion to star formation which can be compared with earlier estimates by 
Kennicutt (1998) and Panuzzo et al. (2003), (see equation 3): 



SFRin 
M yr- 1 . 



4.2 x 10" 44 



erg.s - 



(33) 



This agrees well with these earlier estimates. However, in general the dust shell will be 
more incomplete, and a greater fraction of luminosity is emitted longward of the 100 /im 
band, so that the numerical factor in equation 33 will be larger. From table 4 we find 
that it can range up to 11 x 10~ 44 , so that star formation estimates using equation 3 may 
underestimate the star formation rate by a factor of over two in certain circumstances. 
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Table 4: Model Spizer & IRAS 25/xm luminosities (given in terms of logL/ L & units for a star 
formation rate of l.OM yr _1 ) as a function of both r and P/k. 



Waveband: 


r 


P/k : 








(Myr) 


10 4 


10 6 


10 7 


ZogL MIPS (24) 












1 


7.395 


7.676 


7.952 




2 


7.331 


7.676 


7.976 




4 


7.288 


7.683 


8.000 




8 


7.271 


7.694 


8.022 




16 


7.267 


7.706 


8.039 




32 


7.268 


7.714 


8.049 


logL mAS (25) 












1 


7.739 


8.003 


8.245 




2 


7.682 


8.006 


8.269 




4 


7.645 


8.014 


8.295 




8 


7.631 


8.026 


8.324 




16 


7.630 


8.038 


8.337 




32 


7.630 


8.046 


8.348 
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Table 5: Model IRAS IR and FIR luminosities (given in terms of logL/L Q units for a star 
formation rate of l.OM yr -1 ) as a function of r and P/k. 



Waveband: 


T 

(Myr) 


P/k : 
10 4 


10 6 


10 7 


logL Fm 


1 


9.131 


9.164 


9.063 




2 


9.111 


9.200 


9.155 




4 


9.096 


9.263 


9.261 




8 


9.103 


9.329 


9.342 




16 


9.116 


9.380 


9.450 




32 


9.127 


9.405 


9.496 


logL m 


1 


9.375 


9.444 


9.471 




2 


9.379 


9.525 


9.536 




4 


9.402 


9.556 


9.616 




8 


9.443 


9.626 


9.686 




16 


9.478 


9.681 


9.758 




32 


9.501 


9.707 


9.793 
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11. Conclusions 

In this paper we have attempted to provide a self-consistent theoretical modelling of 
the pan-spectral energy distribution of starbursts within the context of a simplified one 
dimensional evolutionary model of the component H II regions it contains. Although this 
modelling is rather more sophisticated than many models presented hitherto, it still suffers 
from a number of important limitations which we feel should imply a caveat emptor to the 
hopeful observationalist who would like to use the results presented here. Enumerating these 
limitations: 

• The clusters are treated as having a single characteristic mass; 1O 4 M , rather than 
being distributed over all possible masses. 

• The one-dimensional model H II evolution is oversimplified, as it does not take into ac- 
count mixing and dynamical instabilities which certainly will occur in two dimensions, 
nor does it take account of the intrinsic inhomogeneities in the surrounding ISM. These 
effects have been incorporated in a 'correction factor' which reduces the efficiency of 
the stellar winds in inflating the surrounding bubble of ionized gas, and its magnitude 
is chosen so as to provide ionization parameters U in the range observed. Since U also 
determines the dust temperature in the surrounding shell, this choice of the correction 
also ensures that our computed dust temperatures in the molecular shells are more or 
less correct. 

• We appear to be missing a population of compact and ultra-compact H II regions 
around single OB stars or small OB clusters. These will provide an excess of flux in 
the 10 — 30/imband, and would allow a better match to the observed SEDs of starbursts 
and to the observed colorxolor diagrams. 

• The dynamical evolution of the central cluster is not taken into account. This will 
alter the geometry of the stars with respect to the molecular gas, and would lead to 
different dust temperature distributions and different effective covering factors of the 
molecular gas, /, as a function of time. This effect will become more important at 
high P/k because the H II region expansion 'stalls' at an earlier phase of its evolution. 

• Collective effects are not taken into account. That is to say, we have not considered 
mutual overlap of H II regions along the line of sight, such as certainly occurs in very 
dense starbursts such as Arp220. The treatment of such cases would require the SEDs 
computed here to be input as source functions into a complex dust radiative transfer 
code such as used by Popescu et al. (2000). 
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• The diffuse component of the dust has only been included in respect of the attenuation 
of stellar light. The emission from this dust, which is expected to become dominant 
beyond 100/iin (Popescu et al. 2000), is not considered. Furthermore, old stars, that 
is to say stars older than 10 s years old, have not been included. As well as providing 
an additional component of direct light in the visual/NIR spectral range, these stars 
are also important in powering the diffuse dust emission component. 

• Observations of starburst galaxies e.g. Meurer, Heckman & Calzetti (1999); Calzetti 
(2001) show that in the visible and UV, much of the wavelength-dependent attenuation 
by dust can be understood as due to a foreground diffuse screen of material. Much 
of this dust would have to lie outside the main star-forming region. In our models 
we do not take this material into account, although we have provided the expected 
attenuation properties of this material in table 3. This material may well provide an 
important contribution to the far-IR dust continuum beyond 100/xm, and would again 
have to be accounted for within a full galaxian dust radiative transfer code. 

Despite these limitations, these models may well prove useful as providing SED templates 
for individual star formation regions within normal galaxies or for starburst galaxies that do 
not contain an important old star population, or in which there is no extended diffuse dust 
component. For this reason, we are making available electronically both the computed SEDs, 
the attenuation functions and the filter transmissions (sampled at the same frequencies) for 
the various space missions referred to in this paper. These are available in the form of 
tab-delimited spreadsheets in the electronic version of the journal. 

Despite the limitations listed above, this work provides a number of important advances: 

• We have demonstrated how the dust attenuation law of Calzetti (2001) for starburst 
galaxies can arise as a natural consequence of turbulence on photochemistry in the 
diffuse phases of the ISM. In an earlier paper, Fischera et al. (2003) showed how the 
large values of Ry observed in starburst galaxies arise naturally from the log-normal 
column density distribution imposed on the diffuse ISM by turbulent processes. In that 
work, the observed absence of a 2200A dust absorption feature remained unexplained. 
Here we have shown how, if most of the interstellar carbon is locked up into PAHs 
rather than graphitic or amorphous organic grains, then the photodissociation of the 
PAHs in the strong UV radiation fields of starburst galaxies can explain the weakness 
of the 2200A feature in absorption, whilst preserving the strong IR PAH emission 
features originating from the photodissociation regions of molecular clouds. 

• We have demonstrated that two ISM parameters control the form of the SEDs; the 
pressure in the diffuse phase of the ISM (or, equivalently, the density), and the molec- 
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ular cloud dissipation timescale. Long dissipation timescales for the molecular clouds 
enhance the PAH IR emission features by increasing the molecular cloud covering fac- 
tors, but serve to obscure the central clusters in the visible and the UV. High pressures 
put the dense dusty gas closer to the exciting stars giving higher stochastic dust grain 
temperatures, i.epstem We have demonstrated that the models are capable of describ- 
ing the SED of well-observed starbursts over more than three decades of frequency, and 
we have shown that bolometrically-based star formation rate estimates and foreground 
dust attenuations can be easily derived from this fitting. 

• We have shown that the models reproduce the observed IRAS colors and range of 
colors for starburst galaxies, and we have provided diagnostic colorxolor plots for the 
Spitzer Space Observatory MIPS and IRAC instruments which offer the potential to 
solve for the two controlling parameters of the SED. 

• We have provided fitting formulae or tables to enable the conversion of observed fluxes 
to star formation rates in the UV (GALEX), at optical wavelengths (Ha), and in the 
IR (IRAS or the Spitzer Space Observatory) We have shown that the 25/zm fluxes 
are particularly valuable as star formation indicators since they depend little on the 
molecular gas dissipation timescale, and only weakly upon the ISM pressure. 

Notwithstanding the limitations of these models, we hope that these results have been 
presented in such a way as to be useful for observers who are seeking both to understand 
star formation in galaxies, and in the Universe at large. 
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